"Medicine, Faculty of"@en . "DSpace"@en . "UBCV"@en . "Kim, Jae Gak"@en . "2013-07-25T09:09:50Z"@en . "2013"@en . "Master of Science - MSc"@en . "University of British Columbia"@en . "The focus of my investigation is the hyperpolarization-activated cyclic nucleotide-gated cation (HCN) channel, also known as the pacemaker channel. There are four mammalian isoforms (HCN1-HCN4) that share about 60% sequence identity with each other, all activated by hyperpolarization of the membrane potential and permeable to both potassium and sodium. The major differences among four isomers are their responses to binding of cyclic adenosine monophosphate (cAMP), the rate of channel opening and closing, and their dependence on voltage. Recent studies have suggested that the opening and closing of HCN channels involve a step that is voltage independent, which depends upon a region that resides within the S4 and S6 transmembrane region. My study is based upon recent data from the Accili lab in which substitution of a phenylalanine (F) residue near the inner activation gate of the HCN2 channel to alanine (A) dramatically and preferentially slows down channel closing and decreases the dependence of closing on voltage. A 6-state, but not a 4-state, cyclic allosteric model incorporating voltage-dependent transitions moving the channels between resting and active states and voltage-independent transitions between closed and open states was able to describe the complex opening and closing of both the wild type and F/A mutant channels in response to changes in voltage. These models also predicted a significant opening probability between the open and closed states when the channel resides in the resting state."@en . "https://circle.library.ubc.ca/rest/handle/2429/44696?expand=metadata"@en . "A kinetic model for the HCN2 ion channel by Jae Gak Kim B.Sc., The University of British Columbia, 2010 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Cell and Developmental Biology) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) July 2013 \u00C2\u00A9 Jae Gak Kim, 2013 ii Abstract The focus of my investigation is the hyperpolarization-activated cyclic nucleotide-gated cation (HCN) channel, also known as the pacemaker channel. There are four mammalian isoforms (HCN1-HCN4) that share about 60% sequence identity with each other, all activated by hyperpolarization of the membrane potential and permeable to both potassium and sodium. The major differences among four isomers are their responses to binding of cyclic adenosine monophosphate (cAMP), the rate of channel opening and closing, and their dependence on voltage. Recent studies have suggested that the opening and closing of HCN channels involve a step that is voltage independent, which depends upon a region that resides within the S4 and S6 transmembrane region. My study is based upon recent data from the Accili lab in which substitution of a phenylalanine (F) residue near the inner activation gate of the HCN2 channel to alanine (A) dramatically and preferentially slows down channel closing and decreases the dependence of closing on voltage. A 6-state, but not a 4-state, cyclic allosteric model incorporating voltage-dependent transitions moving the channels between resting and active states and voltage-independent transitions between closed and open states was able to describe the complex opening and closing of both the wild type and F/A mutant channels in response to changes in voltage. These models also predicted a significant opening probability between the open and closed states when the channel resides in the resting state. iii Preface Electrophysiological data, protocol, and F431A mutation were provided by Wai Wong, a M.Sc. student in the Accili lab at the University of British Columbia. The construction of 4-state model in ordinary differential equations form and development of analytical solution were done with assistance from Dr. Yue-Xian Li and his Ph.D student Jia Gou from Department of Mathematics at the University of British Columbia. Besides above contribution, all of the work presented is original and unpublished work by the author. The main works include developing the 6-state model, writing MATLAB codes, converting the experimental data, and fitting the 4 and 6-state models to the experimental data. iv Table of contents Abstract ................................................................................................................................................ ii Preface ................................................................................................................................................ iii Table of contents ................................................................................................................................. iv List of tables ....................................................................................................................................... vii List of figures .................................................................................................................................... viii List of equations ................................................................................................................................... x List of abbreviations ............................................................................................................................ xi Acknowledgements ............................................................................................................................ xii Dedication ......................................................................................................................................... xiii 1. Introduction .................................................................................................................................. 1 1.1. The HCN channel and Ih current.................................................................................... 1 1.1.1. HCN channel topology and its regulation ...................................................................... 2 1.1.2. HCN isoforms ............................................................................................................... 5 1.1.3. Opening kinetics and proposed models ......................................................................... 5 1.2. Proposed work ............................................................................................................... 8 1.2.1. Cyclic models ................................................................................................................ 8 1.2.2. F431A mutation: Direct evidence of voltage-independent step ................................... 12 1.3. Rationale and hypothesis ............................................................................................. 15 2. Methods and approaches ............................................................................................................ 16 2.1. Cyclic allosteric models .............................................................................................. 16 2.2. System of ordinary differential equations .................................................................... 17 2.3. Solving the differential equations ................................................................................ 20 2.3.1. Initial conditions .......................................................................................................... 20 2.3.2. Steady-state (analytical) solution ................................................................................. 22 2.3.3. Eigenvalues-eigenvectors solution .............................................................................. 23 2.3.3.1. Eigenvalues and eigenvectors ....................................................................... 24 2.3.3.2. General solution using eigenvalues and eigenvectors.................................... 25 2.4. Conversion of experimental data ................................................................................. 27 v 2.4.1. Conversion of activation protocol ............................................................................... 27 2.4.2. Conversion of deactivation protocol ............................................................................ 30 2.4.3. Conversion of F431A mutant channel data.................................................................. 32 2.5. Fitting of experimental data ......................................................................................... 33 2.5.1. Guess game ................................................................................................................. 33 2.5.2. Genetic Algorithm (GA) of MATLAB ........................................................................ 34 2.5.3. Goodness of the fitting ................................................................................................ 36 3. Results ........................................................................................................................................ 37 3.1. Converted experimental data ....................................................................................... 37 3.1.1. The wild type data ....................................................................................................... 37 3.1.2. F431A mutant HCN2 channel data .............................................................................. 39 3.1.3. Overlapped wild type and F431A HCN2 data ............................................................. 41 3.2. The 4-state model ........................................................................................................ 43 3.2.1. Analytical solution and its reproduced data ................................................................. 43 3.2.2. Wild type HCN2 channel fitting .................................................................................. 47 3.3. The 6-state model ........................................................................................................ 49 3.3.1. Wild type HCN2 data fitting ....................................................................................... 49 3.3.2. 4-state model vs. 6-state model ................................................................................... 51 3.3.3. F431A mutant data fitting ........................................................................................... 55 3.4. Goodness of the fit ...................................................................................................... 60 3.5. Parameters ................................................................................................................... 63 4. Discussion .................................................................................................................................. 65 4.1. Motivation for the mathematical model ....................................................................... 65 4.2. The 4-state model ........................................................................................................ 67 4.2.1. Fitting of channel activation ........................................................................................ 67 4.2.2. Fitting of the deactivation ............................................................................................ 68 4.2.3. Analytical solution ...................................................................................................... 69 4.3. The 6-state model ........................................................................................................ 71 4.3.1. Fitting of the wild type current data............................................................................. 71 4.3.2. The fitting of F431A mutant channel .......................................................................... 76 4.4. Parameters ................................................................................................................... 79 4.4.1. Voltage dependent parameters ..................................................................................... 79 vi 4.4.2. Voltage independent parameters.................................................................................. 81 4.5. Physiological significance ........................................................................................... 85 4.6. Limitations .................................................................................................................. 86 4.7. Future applications and possible improvements .......................................................... 87 5. Conclusion ................................................................................................................................. 88 References .......................................................................................................................................... 89 Appendix ............................................................................................................................................ 92 Appendix A : MATLAB code ........................................................................................................ 92 4-state model simulation with analytical solution ....................................................................... 92 Codes for G.A. ........................................................................................................................... 94 Initial response (customized MATLAB function) ...................................................... 94 dol_myFit (customized MATLAB function) .............................................................. 94 dol_GA (code for setting and running G.A.) .............................................................. 94 6-state simulation for individual analysis ....................................................................................... 95 vii List of tables Table 3.1 Analytical solution of the 4-state model .......................................................44 Table 3.2 Wild type fitting evaluation for the 4-state model ........................................60 Table 3.3 Wild type fitting evaluation for the 6-state model ........................................61 Table 3.4 F431A fitting evaluation for the 6-state model .............................................62 Table 3.5 Parameters for the 4-state model ...................................................................63 Table 3.6 Parameters for the 6-state model (wild type) ................................................64 Table 3.7 Parameters for the 6-state model (F431A) ....................................................64 Table 4.1 Wild type experimental time .........................................................................67 Table 4.2 Percentage changes for the voltage dependent parameters ...........................80 Table 4.3 Percentage changes for the voltage independent parameters ........................81 viii List of figures Figure 1.1 Contribution of HCN channel in the cardiac cell ............................................2 Figure 1.2 Topology of HCN channel ..............................................................................3 Figure 1.3 Effect of cAMP ...............................................................................................4 Figure 1.4 Rate of opening in HCN1 & HCN2 channel ...................................................7 Figure 1.5 Different cyclic models .................................................................................11 Figure 1.6 Location of F431A mutation in HCN2 channel ............................................12 Figure 1.7 Effects of mutation ........................................................................................14 Figure 2.1 The restricted 6-state model ..........................................................................16 Figure 2.2 Overall procedure to get expression for Cr state ...........................................18 Figure 2.3 System of ODEs for the 4-state model ..........................................................18 Figure 2.4 Matrix form of the 4-state model ..................................................................19 Figure 2.5 System of ODEs and matrix form of the 6-state model ................................19 Figure 2.6 Activation (top) and deactivation (bottom) protocols ...................................21 Figure 2.7 Zero in derivatives for the system of ODEs ..................................................23 Figure 2.8 Matrix multiplication ....................................................................................24 Figure 2.9 Finding eigenvalues and eigenvectors ..........................................................24 Figure 2.10 Simple form for the matrix of the 4-state model ...........................................25 Figure 2.11 Activation protocol and its current trace .......................................................28 Figure 2.12 Single current trace for activation .................................................................29 Figure 2.13 Deactivation protocol with current trace .......................................................31 Figure 3.1 Converted wild type data for activation and deactivation .............................38 Figure 3.2 Converted F431A data for activation and deactivation ................................40 Figure 3.3 Overlapped curves from wild type and F431A data .....................................42 Figure 3.4 Reproduced simulations using analytical solution ........................................46 Figure 3.5 Wild type fitting ............................................................................................48 Figure 3.6 Fitting of the 6-state model for the wild type data ........................................50 Figure 3.7 Two possible pathways in opening ...............................................................51 ix Figure 3.8 Changes in individual state at 0mV ..............................................................53 Figure 3.9 Changes in individual state at +20mV ..........................................................54 Figure 3.10 F431A activation fitting with the 6-state model ...........................................56 Figure 3.11 Deactivation fitting with the 6-state model (-20mV) ....................................57 Figure 3.12 Deactivation fitting with the 6-state model (0mV) .......................................58 Figure 3.13 Deactivation fitting with the 6-state model (20mV) .....................................59 Figure 4.1 More hyperpolarizing voltage group activation process for wild type .........74 Figure 4.2 Less hyperpolarizing voltage group activation process for wild type ..........75 Figure 4.3 More hyperpolarizing voltage group activation process for F431A .............77 Figure 4.4 Less hyperpolarizing voltage group ..............................................................78 Figure 4.5 Modified parameters .....................................................................................79 Figure 4.6 Deactivation simulation of two models ........................................................83 x List of equations Equation 2.1 General solution in terms eigenvalues and eigenvectors ..............................26 Equation 2.2 Formula for fitting evaluation .......................................................................36 xi List of abbreviations A: Alanine Ca: Active closed state cAMP: Cyclic adenosine monophosphate CNBD: Cyclic nucleotide binding domain Cr: Resting closed state D.N.A.: Deoxyribose nucleic acid F: Phenylalanine F\u00E2\u0088\u009E: Infinite factor G.A.: Genetic Algorithm MATLAB: Matrix Laboratory MC: Max-current LC: Leak-current Oa: Active open state ODE: Ordinary differential equation Or: Resting open state HCN: Hyperpolarization-activated cyclic nucleotide-gated cation H-H: Huxley-Hodgkin If : Funny current Ih : Hyperplarization-activated current SpIH: Sea urchin HCN xii Acknowledgements Thank you for supporting my research in every possible way from my lab. A big thank-you goes to Dr. Eric Accili for his endless motivation and academic inspiration to my research. He is one I truly respect and admire. Thank you, Mr. Leo Ng, for his countless corrections and support as he refined English usage and biochemical facts in this thesis. The electrophysiological data and help on software & computer issues are credited to Mr. Wai Wong. Thank you to the rest of my lab members as well. Thank you, Dr. Yue-Xian Li and Ms Jia Gou for keen mathematical consultation and discussion. Without their help, I would not have started and finished this research. Thank you to the faculties and staff of the CELL program who make it possible for me to study and research at UBC. Thank you to my family and friends who supported me along the way, physically, mentally, and financially. xiii Dedication To. Mina & Mom 1 1. Introduction 1.1. The HCN channel and Ih current As its name implies, HCN channel is voltage-gated ion channel that is responsible for producing a current called namely If or Ih. Both names, \"funny current\" (If) and \"hyperplarization-activated current\" (Ih), are in reference to the properties of the current that, upon their discovery, were found to be different from many of the other known current systems (Accili et al., 2002; Baruscotti et al., 2005; Biel et al., 2009). Unlike many voltage- gated ion channels that are activated and open when they are depolarized, HCN channels are slowly activated open by hyperpolarization of the membrane potential. Upon opening, both Na+ and K+ ions conduct, although the latter is more permeable than the former. The general role for the current is to depolarize the membrane potential following hyperpolarization and, thus, its activation often leads to the threshold level for the generation of the action potential in many different excitable cells. Therefore, the physiological role of HCN channels is not only limited to pacemaking activities, but also extended to neuronal excitability, dendritic integration, and network oscillations (Biel et al., 2009; Baruscotti et al., 2010; Benarroch, 2013). Figure 1.1 shows the contributions from different ion channels to the action potential for pacemaker cells of the sinoatrial node. The sinoatrial HCN channel opens upon hyperpolarization of the action potential beyond ~ -40mV, voltages at which the net movement of cations is inward. Although it is not the only contributor (Lakatta & DiFrancesco, 2009), the opening sinoatrial HCN channels impacts both the diastolic depolarization phase of the action potential and heart rate (DiFrancesco, 1993). 2 Figure 1.1 Contribution of HCN channel in the cardiac cell: Contribution of different ion channels to generate action potential at sinoatrial nodal cell. HCN channel is activated around -55mV and initiates the depolarization for threshold level of action potential. Source: (Sherwood & Cengage Learning (Firm), 2010) 1.1.1. HCN channel topology and its regulation The opening of HCN channels is thought to physically alter the structural components of the channels (Biel et al., 2009) and, therefore, an understanding channel structure and topology is important for understanding the HCN gating mechanism. The proposed topology of HCN channel is based on that of potassium channel, another member of the voltage-gated ion channel (Accili et al., 2002). As shown in Figure 1.2, HCN channels are tetramers and each of the 4 subunits is composed of three domains: an N-terminal domain, 6 transmembrane segments core (S1-S6), and a C-terminal domain. Although S5 and S6 form the central pore of the channel, the S4 and C-terminus are important in control of channel opening and closing. 3 In HCN channels, it is thought that the positively charged voltage sensor, S4, detects hyperpolarizing potential and moves towards the intracellular side of the cell, the same way as that for depolarization-activated potassium channels, but this movement is inversely coupled to the opening of the pore (Mannikko et al., 2002; Bell et al., 2004; Vemana et al., 2004). Furthermore, it has recently been suggested that the coupling of voltage-sensor to the pore is weak and that it may act more as a modulator of ion flow (Ryu & Yellen, 2012). Figure 1.2 Topology of HCN channel: Proposed topology of HCN channel being shown as tetramer and its subunit composition: transmembrane domains (S1-S6), voltage sensor (S4), and cyclic nucleotide binding domain (CNBD). Source:(Postea & Biel, 2011) 4 The C-terminus contains the cyclic nucleotide binding domain (CNBD), where cAMP directly binds and enhances channel opening. The binding of cAMP shifts the range over which the HCN pore opens to less negative potentials (DiFrancesco & Tortora, 1991). Figure 1.3 shows the effect of saturating amount of cAMP (100\u00C2\u00B5M) on channel opening probability, and compares to that in the absence of cAMP; this experiment is carried out in a patch excised from an oocyte expressing the HCN2 isoform. In this isoform, and under these conditions, the necessary potential to open 50% of a given population of total HCN channels (V1/2) is -125\u00C2\u00B12.0mV in control. However, 100\u00C2\u00B5M of cAMP increases the rate of activation, and V1/2 becomes-111\u00C2\u00B11.7mV (Craven & Zagotta, 2004). In the HCN2 isoform, cAMP also increases the maximum amount of current, a result that may reflect an increase in channel number or maximum open probability. The fundamental mechanism of this cAMP- dependent modulation of HCN channels is \u00E2\u0080\u009Cdeinhibition\u00E2\u0080\u009D. CNBD, behaves as an auto- inhibitory domain to lower the opening probability, and the addition of cAMP removes this tonic inhibition of CNBD and increases channel activity (Wainger et al., 2001). Figure 1.3 Effect of cAMP: Activation curve of HCN2 channel with (red) and without (black) cAMP. Normalized conductance in y and voltage in x-axis. Source:(Craven & Zagotta, 2004) 5 1.1.2. HCN isoforms There are four known mammalian isoforms, HCN 1-4 (Ludwig et al., 1998; Santoro et al., 1998; Ishii et al., 1999). All of them conduct the Ih and have similar permeation properties, but there are major differences to binding of cAMP and rate of activation (Biel et al., 2009). The HCN 1 channel opens and closes the fastest, followed by HCN2, HCN3 and HCN4, in that order. The HCN1 channel is not greatly affected by cAMP binding and the HCN3 channel is not affected at all by binding of this molecule. HCN 2 and 4 have about same positive shift from binding of cAMP, the primary channel to be studied in this research is HCN2 channel. 1.1.3. Opening kinetics and proposed models The opening process of HCN channel, which is activated by hyperpolarization and enhanced by binding of cAMP, is believed to involve multiple complex steps according to experimental observations and data (DiFrancesco, 1999; Altomare et al., 2001; Mannikko et al., 2005; Craven & Zagotta, 2006). DiFrancesco observed that the current traces of channel opening have a distinctive initial delay, which gives a sigmoidal shape in the current trace (DiFrancesco, 1984). This specific feature is difficult to reconcile with the Hodgkin-Huxley (H-H) model even though the H-H model has successfully described the kinetics of most membrane currents from voltage-gated ion channels in the heart and other excitable tissues since 1952. He also demonstrated that initial delay can be described by the third power of exponential function, but post-delay current trace of the opening process can be fitted by a single power exponential function. These findings led to the conclusion that the opening process must be constructed with different steps, and explained why no single mathematical 6 expression can describe the current trace of channel opening process simultaneously. This requires the involvement of a kinetic model, which is composed of multiple mathematical equations and each of them may be accountable to each of the different steps in channel opening process. Currently, there is experimental evidence to support the fact that the channel opening process may be governed by two different steps: voltage-dependent and voltage-independent steps (Chen et al., 2007). The involvement of voltage-dependent step in the opening process is most likely required since the HCN is a voltage-gated ion channel that is activated upon hyperpolarization. For the existence of voltage-independent step, the evidence can be found from Figure 1.3. The open probability curve forms the plateau shape in the absence of cAMP at extreme negative potential around -150mV. This suggests that negative potentials beyond -150mV do not open more HCN channels. In other words, the voltage-dependent step of channel opening process is saturated and no significant change is induced even with more negative potentials. Nevertheless, exposure to 100\u00C2\u00B5M of cAMP increases the maximum conductance, suggesting that cAMP enhances another step that does not depend upon voltage. Another evidence for the presence of a voltage-independent step can be found in the Figure 1.4, which contains the rates of opening for HCN1 and HCN2 channels. As shown, at a voltage more negative than -170mV, the rate of opening changes little despite the applied potentials for HCN2 channel. However, HCN1 channel does not exhibit such a voltage- independent behavior. To further investigate on this voltage-independent phenomenon, S4- S6 domains of HCN1 were substituted with that of HCN2 channel. HCN1 channel with S4- S6 domains from HCN2 channel also exhibits the same behaviour. Original authors 7 concluded that S4-S6 domain must be responsible for the voltage-independent step (Chen et al., 2007), but a more specific location of the voltage-independent step within the channel has not been identified. Finally, there is evidence to suggest that there is voltage-independent current carried by HCN channels (Gauss et al., 1998; Ishii et al., 1999; Chen et al., 2001; Proenza et al., 2002; Macri & Accili, 2004). For HCN2 and the sea urchin HCN channel SpIH, this voltage- independent current is estimated to be around 2% and 8% of the fully-activated current, respectively (Proenza & Yellen, 2006). Support for current flow through the resting channel has been recently provided by measurement of gating current in locked closed and locked open SpIH channels, in which an open probability of ~17% was estimated at rest (Ryu & Yellen, 2012). Figure 1.4 Rate of opening in HCN1 & HCN2 channel: The rate of opening (1/tau in y-axis) between HCN1 and HCN2 at different potentials (mV in x-axis). HCN2 shows less voltage sensitivity in the rate after -160mV but HCN1 does not have such a pattern. Source: (Chen et al., 2007) 8 1.2. Proposed work The proposed work is mainly to improve and construct the kinetic model to successfully describe both voltage-dependent and voltage-independent steps in HCN channel opening from previously published models. The simplest model is the sequential model but it does not completely reproduce channel behaviour because it does not allow channels to open when they are unbound or partially bound with cAMP; HCN channels can be readily opened without cAMP. The next progression is the Monod, Wyman, and Changeux (MWC) model, which is originally developed to characterize the allosteric behaviour of hemoglobin. MWC model well describes the dual modulation by hyperpolarization and binding of cAMP in which a concerted conformational transition occurs from the closed to the open state of the channel with any level of extracellular cAMP (DiFrancesco, 1999). However, the model is limited by the assumptions that there are four identical binding sites (Mannikko et al., 2005; Craven & Zagotta, 2006). Moreover, MWC model does not necessarily involve the voltage- independent step that needs to be proved. This leads to the cyclic model that contains both voltage-dependent and independent steps. 1.2.1. Cyclic models The present research was inspired by Siegelbaum\u00E2\u0080\u0099s 4-state cyclic model, which was used to explain the kinetic behaviour of channel opening and closing in HCN1 and HCN2, and focused on voltage-independent behaviour (Chen et al., 2007). The model consists four states, connected by voltage-dependent and voltage-independent steps (Fig 1.5). The four states are divided into open (O) and closed (C) states, denoted as resting (Or, Cr) and active (Ca,Oa) sub-states. The horizontal transition between states is voltage-dependent step, where 9 that the voltage sensor movement in response to changes in membrane potential is predicted. The respective movement of voltage sensor switches resting state to active sub-states. The vertical transitions between closed and open, is a voltage-independent step. In Chen et al, a linear 3-state model was first attempted before the cyclic 4-state model. However, the simpler linear model was not able to reproduce the experimental data for channel closing. The next natural approach is simply adding one more state. A cyclic 4-state model provided a better fitting result in both opening and closing. Beside the 4-state model, there are different cyclic models that explain and represent HCN channels at a molecular level. From the topology of the HCN channel, there are four subunits and each subunit has its own voltage sensor. The simple 4-state model assumes all four voltage sensors move at same time in a single voltage-dependent step. However, several studies have utilized models with additional states that account for movement of each voltage sensor in every subunit. The next simplest model after 4-state is the 6-state model, also referred to as a dimer of dimers model (Craven & Zagotta, 2006; Bruening-Wright & Larsson, 2007). It has two voltage-dependent steps in opening or closing process and only two of the four voltage sensors move in single voltage- dependent step. The most detailed and complete model is the 10-state model in which individual voltage sensor movement is illustrated (Altomare et al., 2001; Ryu & Yellen, 2012). The more complex the model is, the better the model will be in describing the features of data in more details and greater accuracy. However, the purpose of my research focuses on the simplest model that can be used for fitting and for fully understanding the role of each parameters and states in the model as much as possible. Therefore, I began with the 4-state model but found that a more complex model (Fig 1.5) was necessary to properly show 10 typical characteristics of HCN channel gating, especially for the sigmoidicity in the beginning of both opening and closing. 11 Figure 1.5 Different cyclic models: 3 different cyclic models are listed. A and B are mainly used in research. 10-state model (C) is provided to understand previous attempt to describe gating of HCN channel. A) 4-state model B) 6-state model (Dimer of Dimers) C) 10-state model and its molecular representation Source: (DiFrancesco, 1999) 12 1.2.2. F431A mutation: Direct evidence of voltage-independent step Originally, F378A mutant was created in the HCN1 and HCN2 isoforms from our lab, as part of an alanine screening of the inner pore. Figure 1.6 displays the exact location of the mutation, at the pore within the S6 segment. Figure 1.6 Location of F431A mutation in HCN2 channel: Exact location of point mutation at S6 region of HCN2 channel. This single mutation changes amino acid from F to A. Source: (Chan et al., 2009) In both mutants, a dramatic delay in closing time was observed, which did not appear to depend strongly on voltage. Interestingly, the activation curves of the HCN2 mutant and wild type channel are not largely different despite of the fact that deactivation of wild type and mutant channel is very different. Together with similarity in rate of opening and the similarity in activation curves (Fig 1.7) suggest that channel opening is not as 13 influenced by the mutation as is closing. The experimental data on the HCN2 channel expressed in the oocyte and measured using the two microelectrode approach are shown in Figure 1.7. It is true that closing time in both wild type and F431A mutant channel have voltage independent behaviour as the applied voltage becomes greater than 0mV. However, the important message is that F431A mutant channel has such a plateau pattern starting at - 40mV and maintains it throughout the whole voltage range during deactivation. This enhanced voltage independent behaviour is easier to visualize when the experimental data is converted to open probability. Therefore, F431A mutant is a more direct evidence to prove the presence of voltage-independent step along with mathematical model that correctly characterizes it. 14 Figure 1.7 Effects of mutation: Effect of F431A mutation in activation curve (A) and rate of opening (B). No significant shift is observed in activation curve but closing of mutant channel become very delayed and less sensitive to the voltage. A) Effect in activation curve B) Effect in rate of opening Source: Wai Wong, MSc. Student in Accili Lab 15 1.3. Rationale and hypothesis Based on the background discussed above, it seems reasonable to propose that the F to A mutation in the HCN2 channel disrupts a voltage-independent transition during channel closing. I hypothesize that the 4-state or 6-state model with appropriately modified parameters will successfully describe both wild type and mutant HCN 2 channel data. Furthermore, this finding will be used as validation and direct evidence for the presence of a voltage independent step at its responsible location (inner pore region at S6). 16 2. Methods and approaches 2.1. Cyclic allosteric models I utilize a 4-state model and the restricted 6-state model, to simulate the currents of the HCN2 channel. The 4-state model is the same as shown earlier in Figure 1.5; I first validated this model and then used it to simulate the experimental data from our own experiments (from Wai Wong). I also utilized a 6-state model, keeping in mind that one of the primary goals is to build the simplest mathematical model that would successfully explain HCN channel data. Figure 2.1 shows the restricted 6-state model, in which the transitions between Cr and C1, and that between C1 and C2 are identical. The transitions between O2 and O1, and that between O1 and Or, are also the same. From here on after, the 6-state model will be referred as the restricted 6-state model. Figure 2.1 The restricted 6-state model: Diagram of restricted 6-state model and direction of given parameters. Horizontal transition is voltage-dependent and vertical transition is voltage-independent step. \u00CE\u00B1, \u00CE\u00B1\u00E2\u0080\u0099 \u00CE\u00B2, \u00CE\u00B2\u00E2\u0080\u0099 g, g\u00E2\u0080\u0099,g'' h, h\u00E2\u0080\u0099,h'' Direction \u00E2\u0086\u0092 \u00E2\u0086\u0090 \u00E2\u0086\u0093 \u00E2\u0086\u0091 17 2.2. System of ordinary differential equations The first task that must be accomplished was converting the 4-state model into mathematical language. Since the model itself is supposed to represent the population or probability of ion channels at each state at any given time, usage of the differential equation is the most appropriate. The differential equation computes the output of the system by using time as the input. The developed differential equations, which mathematically reproduce the 4-sate model, only have time as their variable. It is true that 4-state model has voltage parameters but those parameters will not be considered as input variables like time. These voltage parameters are dependent on voltage and because voltage changes, it can be viewed as another input variable. However, actual opening and closing processes were performed at fixed voltage values so voltage dependent parameters should be treated as the constant terms. Therefore, the proposed differential equations of 4-states model will be a single variable function, more specifically, a first-order ordinary differential equation (ODE). ODE will describe output variables, namely Cr, Ca, Oa, and Or, at any given time like 4-states model. It is relatively easy to configure such a system of ODEs as one carefully follows the model. I will begin with state Cr as the initial starting point. It is clear that Cr is connected to the Ca and Or states and there are a total of four parameters involved in the connection. They are \u00CE\u00B1, \u00CE\u00B2, g\u00E2\u0080\u0099, and h\u00E2\u0080\u0099. Whether the parameters are denoted with prime or not, \u00CE\u00B1 are responsible for transition from left to right and g is for the downward transition. In a similar manner, \u00CE\u00B2 is for transition from right to left and h is for the upward transition. Since there are four parameters with two states, population of Cr is constantly changing and it can be either positive or negative according to the size of parameters and population in each of the connected states. 18 The Cr state redirects its population to Ca and Or states in the magnitudes of \u00CE\u00B1 and g\u00E2\u0080\u0099 respectively Cr). At same time, Cr gains its population from Ca and Or states at the magnitude of \u00CE\u00B2 and h\u00E2\u0080\u0099 respectively ( Ca + Or). The total and net change of Cr state is simply the addition of these expressions of gaining and losing. This new combined net expression explicitly describes the total net change of Cr state in terms of the connected states and parameters. Figure 2.2 summarizes the derivation of the expression for Cr state; a similar approach was used to construct the other states. Figure 2.3 is the final form of the 4- state model in system of ODEs. Figure 2.2 Overall procedures to get expression for Cr state Ca + Or Ca + Or Figure 2.3 System of ODEs for the 4-state model \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B1 \u00CE\u00B2 In order to use the MATALB to solve the system of ODEs, equivalent matrix form is required. Figure 2.4 is the equivalent matrix form of the 4-state model. 19 Figure 2.4 Matrix form of the 4-state model \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B1 For this research, the 4-state model is not the only one mathematical model that has been used. Since the 6-state model is used, system of ODEs and matrix form of it is also required. Again, the derivation is the same as for the 4-state model. Figure 2.5 has both system of ODEs and matrix form of the 6-state model. Figure 2.5 System of ODEs and matrix form of the 6-state model \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B1 20 2.3. Solving the differential equations There are two main generic ways to solve the differential equations. First is finding numerical solution using numerical methods. Second is finding general solution. From those two main methods, there are many different approaches one can take to solve the differential equations. Regardless of the approach, solution should be and would be identical at the end. However, in order to properly and effectively solve any differential equation, appropriate initial conditions need to be assigned. Validation of initial conditions is essential when solving or estimating many mathematical solutions. If one cannot provide acceptable initial conditions, the solution will most likely be wrong, regardless of methods. Therefore, the primary task is to acquire reasonable initial conditions for each of four states in HCN2 channels. 2.3.1. Initial conditions It is important to define the meaning of initial conditions for HCN2 channels. There are two main voltage protocols that measure the flow of ionic current in HCN2 channels: opening (activation) and closing (deactivation) protocols as shown in the Figure 2.6 21 Figure 2.6 Activation (top) and deactivation (bottom) protocols A) Activation protocol: Typical voltage protocol for activation. X-axis is time in seconds and y-axis is the holding potential in mV. Initially starting at very positive voltage (0mV) and change to negative voltage (ex: -110mV) to activate channel to record current. B) Deactivation protocol: Typical voltage protocol for deactivation. X-axis is time in seconds and y-axis is the holding potential in mV. Initially starting at very positive voltage (0mV) and change to negative voltage (ex: -100mV) to activate channel. Then it goes back to more positive voltage (ex: +20mV) to observe current changing during deactivation. For activation, voltage is initially held at a very depolarizing voltage (for example, 0 mV) such that it is assumed most channels are closed. After that, a hyperpolarizing voltage is applied to activate and open the channels. For the case of activation, the initial conditions are individual population of the four states (Cr,Ca,Oa,Or) initially at the holding potential of 22 0mV. On the other hand, fully activated channels are attained in the closing protocol by a hyperpolarizing voltage (for example, -120mV) at the beginning, and these channels are subsequently closed via the application of depolarizing voltages. For the case of deactivation, the initial conditions are the population of each of the four states at -120mV. It is certainly challenging and nearly impossible to collect any measurable ion channel recording to represent the population of the each four states for the construction of the model. These states are likely theoretically proposed characteristics to explain the typical behavior in experimental data of HCN channels. For such a technical limitation, initial conditions have to be calculated or estimated with valid assumptions. In the case of the 4-state model, initial conditions have been calculated using analytical solutions. For the 6-state model, initial conditions are estimated without calculations. However, initial conditions of this specific system are two extremes as the population should be mostly open at very negative voltage, -120mV, and mostly closed at more positive voltage, 0mV. Therefore, estimated initial conditions do not harm the overall quality of solution in a considerable manner for the 6-state model even though they are not explicitly calculated. 2.3.2. Steady-state (analytical) solution One of the most important aspects and characteristics of an initial condition is that it stays as constant value. The populations of the four states should not fluctuate over the course of times at the initial voltage level as they are supposed to hold for long time to ensure that the steady-state of the channel is achieved. From the perspective of the system of ODEs, reaching steady-state means all the derivatives in the equations will become zero. This new change to the system of ODEs is reflected in the Figure 2.7 23 Figure 2.7 Zero in derivatives for the system of ODEs \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B1 \u00CE\u00B2 The derivative, which is taken with respect to time, serves as the rate of population change in each state. As the derivative becomes zero, rate of change in each state also becomes zero at steady-state. System of ODEs can be written as Figure 2.7 and now each state can be solved and express in terms of the other states and parameters because the system is linear. This final expression is called an analytical solution and it determines the exact population in each of the four states during steady-state condition at any holding voltage values. 2.3.3. Eigenvalues-eigenvectors solution Since the initial conditions of each state can be theoretically calculated using steady- state solution, it is now possible to solve the system of ODEs. As previously stated, there are two ways to solve differential equations: numerical solution and general solution. For the purpose of our aims and goals, it is necessary to study meaning of different parameters and states. For such reason, general solution needs to be solved in terms of eigenvalues and eigenvectors. In order to solve and understand general solution, the concept of eigenvalues and eigenvectors is required. 24 2.3.3.1. Eigenvalues and eigenvectors Let\u00E2\u0080\u0099s assume there is an arbitrary square matrix A multiplied with another arbitrary matrix x. A new matrix is produced as matrix B. With specific entries of matrix A, it behaves like a linear transformer on vector x and the effect of matrix A is reproduced as matrix B. For example, if A is and x is , then B will be as shown in the Figure 2.8 Figure 2.8 Matrix multiplication As Figure 2.8 shows, matrix A makes matrix x two times bigger. It is possible to expand and reduce matrix B. When the sign is reversed, it is interpreted that the direction of the matrix is reversed as well. The matrix A behaves exactly like linear transformation of the function. We denoted the direction component as eigenvector and the magnitude component as eigenvalues. It is usually mathematically expressed as x, where x is the eigenvector and is the eigenvalue of A. Figure 2.9 will show the step-by-step derivation to finding eigenvalues and eigenvectors for a given vector as it is defined in x. Figure 2.9 Finding eigenvalues and eigenvectors x 0 0, where I is the identity matrix. In order to find the eigenvalues ( , one needs to find the determinant of when it becomes zero. Once the specific value of , satisfying the determinant of 25 , has been calculated, eigenvectors can be computed after substituting the values back into the original equation. It is very long and complicated even for a 4 x 4 matrix of A, in the case of 4-state model, if it were computed out manually. For this reason, built-in commands of MATLAB are used to calculate eigenvalues and eigenvector. 2.3.3.2. General solution using eigenvalues and eigenvectors For the sake of simplicity, new and simpler notations are used to describe matrices of the 4-state model in Figure 2.10 Figure 2.10 Simple form for the matrix of the 4-state model Now, finding the solution of matrix becomes finding a unknown function that satisfies the equation Generally A will be n x n matrix and in this case of 4- state model it will be a 4 by 4 matrix. For 6-state model, the A will be a 6 by 6 matrix. The 26 general solution of the differential equations from matrix forms can be written as equation 2.1 using eigenvalues and eigenvectors, Equation 2.1 General solution in terms eigenvalues and eigenvectors Equation 2.1 is the exact solution for a system of ODEs in the matrix form, in which terms are constants, terms are eigenvectors, and terms are eigenvalues. terms are calculated using initial conditions. For the 4-state model, one needs to set n=4 and expand the above equations until n=4. Note that is actually 4 by 1 vector array so each of the 4 solutions in represents each of the 4 states. If n=6, which reflects 6-state model, becomes 6 by 1 vector array. 27 2.4. Conversion of experimental data The successful general solution provides the activation and deactivation simulations with modified parameters at every voltage step. One can simply observe and predict the dynamics in the population of each state as voltage is altered. However, the produced graphs are simply solutions of some mathematical equations and they are biologically meaningless unless they correspond to experimental data. In order to extract any relevant information from the simulation, one needs to fit and compare the solution to the actual experimental data. Our experimental data from electrophysiology recording produces current traces as its final product. This means the current traces have to be converted into normalized open probabilities, which are the immediate outputs of the simulation. Fitting to experimental data is mandatory to provide physiological meaning to the simulation. Furthermore, the conversion of experimental data in currents to the normalized open probability is needed for the fitting and drawing physiological conclusions related to mathematical models. 2.4.1. Conversion of activation protocol As previously discussed, there are two protocols and the conversion of activation data is explained first. I begin with the activation protocol along with currents in Figure 2.11. In this Figure, only one particular voltage value and its corresponding current trace is highlighted in the red. However, every other voltage values and their current traces have been converted. The most hyperpolarizing voltage used is -120mV up to -50mV in increments of 10mV. Thus, a given set of data should have a total of 8 useable current traces that can be converted to open probability. Figure 2.12 is the single current trace for -100mV from actual experimental data. 28 Figure 2.11 Activation protocol and its current trace: This Figure is composed of two graphs top one is current trace and bottom one is voltage protocol used to get current reading. X-axis is time in second. Top y-axis is current in uA and bottom y-axis is voltage in mV. 29 Figure 2.12 Single current trace for activation: Single current trace that is obtained from 0mV to -90mV and then back to +35mV. Y-axis is current in uA and x-axis is time in second. Important part of trace is labelled for calculating open probability. The basic equation to convert current to the open probability is \u00E2\u0088\u009E . It is crucial to note that even though leak current is subtracted, it was reported that HCN channels are never closed; channels are always open 2-5%. I arbitrarily set up for the open probability as 2% at fully closed state. The normalized factor, F, will be important as it determines the overall open probability accordingly to the voltage values. For example and as the example Time (second) C u r r e n t (u A ) Infinite factor (F\u00E2\u0088\u009E) Ih Max-current (MC) Leak-current (LC) 30 only, -120mV will have a factor of 1 but -80mV will have factor of around 0.5 if V1/2 is around -80mV. Out of 5 sets of experimental data, 3 of them are used and averaged. The reason for not using all the data is because some data have notable noise and unstable holding potential. 2.4.2. Conversion of deactivation protocol The conversion of the deactivation protocol has is similar to conversion of the activation data. The main differences are that only tail currents are considered and that the closing current traces will be normalized to themselves. In other words, the converted open probability starts from 1 when fully opened, and goes to 0.02, when fully closed regardless of the applied voltages. This way, each open probability curve has the same initial and final points but their rates should be different since closing of the channels is voltage specific. Although three voltage values used for deactivation, -20mV, 0mV, and +20mV, I did not use current traces at -20mV because they were very small and variable in size and rate of change. Thus, only current trace from 0mV and +20mV are converted. Figure 2.13 is a typical deactivation protocol along with its current traces and a zoomed-in representation of the tail current section for the conversion calculation. 31 Figure 2.13 Deactivation protocol with current trace: Top part is deactivation protocol and its recording current traces. Bottom part is one of tail currents being zoomed- out. Probability conversion of deactivation only uses tail current and it has been normalized to itself. Time (second) C u rr en t (u A ) Current forced to be 2% Time = 0 as it is start of conversion Maximum point where open probability is 1 32 Converting deactivation currents into open probabilities is simpler. Each current is normalized to the its own max current, where time is 0 and the smallest current is fixed at 2%, the same arbitrary value for the activation protocol. There is no voltage specific factor such as the infinite factor for the activation. Therefore, each converted open probability of different voltage values starts at 1 (100%) and decreases until it reaches 0.02 (2%) at the end. The assumption of 2% is only applicable if and only if the end of tail current reaches a steady-state. Although every different voltage value has same initial and final values, their time kinetics should reflect the voltage sensitivity. For example, the rate at which the population approaches 0.02 (2%) is faster for a more positive voltage. 2.4.3. Conversion of F431A mutant channel data The conversion process for F431A data is exactly the same for both activation and deactivation. However, because closing time for F431A mutant channels is very slow, the tail currents do not seem to reach steady-state. It is difficult to determine the channel open probability at the end of the experiment. The end values for those cases are assigned with the aid of a simulation.. Instead of 2%, -20mV ends at 9%, 0mV 8%, and +20mV 7%. 33 2.5. Fitting of experimental data The fitting of wild type data is a requirement before fitting with that of the F431A mutant channel. If the mathematical model does not reasonably fit the wild type data, there is no connection and meaning in comparing the phenotypical changes between the mutation and wild type. Unfortunately, the original parameters used in previous model (Chen et al., 2007) do not fit the data from our lab; this may be due to the different experimental procedures for the two labs are different. For this reason, it is necessary to find the right parameters by fitting the wild type data first. Parameters can then be modified to fit the mutation data. Two different methods were used for the fitting. 2.5.1. Guess game First method that was used extensively is called \u00E2\u0080\u009CGuess and Run\u00E2\u0080\u009D. This method seems oversimplified but it is the most natural method in providing a user with the details of each individual parameter. One can view the fitting as shaping and sculpting the curve into a desirable form. One might want to have a smoother curve or faster dynamics over the course of simulation time or any characteristics of the curve in order to reproduce the target data. Each of the eight parameters is responsible for a certain part of the curve. Thus, in order to fit the data more accurately, the specific region modified by a respective parameter must be mapped. By altering the values of these parameters and running the simulation, the after- effect on the curve can be observed. To begin this method, countless random changes are applied to fit the experimental data. I had little success in determining the engaged parameters for specific adjustments within the curves. I therefore expended a considerable amount of effort and time to eliminate 34 parameters that did not greatly impact the fitting. It is true that all eight parameters contribute to the simulation but their magnitudes of contribution are different. The fundamental idea is that if the number of candidate transitions is reduced, the remaining parameters can be solved more easily. Moreover, the ultimate goal is spotting the parameters that induce the linked results in eigenvalues and eigenvectors, since each simulated curve is actually a solution composed of eigenvalues and eigenvectors in a system of ODEs. Several different MATLAB codes have been written to test the degree of alteration in eigenvalues and eigenvectors from changing the parameters. 2.5.2. Genetic Algorithm (GA) of MATLAB Genetic Algorithm (G.A.) is a built-in function from the Optimization Tool Box of MATLAB software. The G.A. tool mathematically calculates the best fitted individual (solution) and discards the less fitting individuals from the pool. Users can set up specific fitness functions that determine the environment solutions have to adapt to. They are also able to set up the possible maximum number of populations and generation. They can even insert mutations or cross-over, if they wish to do so. For any kinetic models that are used for fitting, the G.A. finds the best parameters from 0 to positive infinity on any given curve. Here, curves are converted to open probability from each voltage value. The particular code developed for research does not support multi- fitting. This means fitting of data from different voltage values is done separately to determine the eight parameters in the case of the 4-state model. However, from previous fitting experience and knowledge of modelling, the suitable range of values for voltage- independent parameters are estimated. At the end, only four voltage-dependent parameters 35 remain to be fitted for every single voltage value. Since the voltage-dependent parameters are in exponential form ( , where c and d are constant values) the ultimate goal would be to find each constant terms ( ) and exponent terms ( ). Another important consideration is that voltage-dependent parameters only increase or decrease in one direction. For example, only decreases as voltage increases. However, the best-fitting parameters do not have such a trend and values incorporated into the program for the sake of the best-fitting. The program does not realize and does not care about such a restriction. The only goal of the program is to find the best values. Such a parameter set does fit almost perfectly only for individual voltage. Since developed code does not support a multi-fitting function, multiple simulations are run repeatedly at a single voltage after manually selecting for the most reasonable values of parameters from the code. Once all of values are picked from -120mV to -50mV, re-fit them into exponential form for the single voltage-dependent parameter. The same procedure is repeated for another voltage-dependent parameter until every voltage-parameter has an exponential expression. For the 6-state model, same procedure is used. All voltage- independent parameters are estimated first and then only four voltage-dependent parameters are fitted individually. Fitness function that has been used for the simulation is coded to calculate the differences between simulation and experimental data. The function will consider such a difference as \u00E2\u0080\u009Cfitness\u00E2\u0080\u009D of one specific generation. If G.A. thinks value of fitness is acceptable, it will evaluate that generation to next level and repeat the process until desirable level of fitness is achieved. The fitness function also includes the constraint of total probability being 1 (Total probability = Sum of all open and closing states). If this condition is violated, fitness 36 function will give a huge penalty in fitness for that particular generation. Thus, penalized generation will have poor fitness value so it will be eliminated from the evolution. 2.5.3. Goodness of the fitting It is very difficult to tell the level of fitting by simply looking at it. To quantify the goodness of fit, equation 2.2 was used. The formula is root-mean-square deviation, which is frequently used to measure the difference between prediction and measurements. The smaller the difference, the better the fit is. Equation 2.2 Formula for fitting evaluation , where is the experimental value and is the reproduced value. 37 3. Results A number of graphs in the results show open probabilities. Every graph has an unitless y-axis depicting open probability and a x-axis for time. Most x-axes are in , meaning every 5000 units equals 1 second in real time. The only exception is Siegelbaum\u00E2\u0080\u0099s data (Chen et al., 2007), which is reproduced where the x-axis is in milliseconds. 3.1. Converted experimental data The converted open probability from both the wild type and the F431A mutant data are presented and compared. 3.1.1. The wild type data Figure 3.1 is the converted open probability from the wild type data for both activation and deactivation. As expected, there is faster kinetics for the opening process at more hyperpolarizing voltages and faster closing process at depolarizing voltages. 38 Figure 3.1 Converted wild type HCN2 data for activation and deactivation: Converted open probability during activation and deactivation for wild type HCN2 channel. Holding potential is individually labelled. As expected more negative potentials have faster opening kinetic during activation and more positive potentials have faster closing kinetic during deactivation. 39 3.1.2. F431A mutant HCN2 channel data Figure 3.2 is the converted open probability from the F431A mutant HCN2 channel data for both activation and deactivation. As it is expected, the open probability curves for activation is not significantly different from that of wild type. However, open probability curves for deactivation definitely show delayed and stronger and more prominent voltage independent behaviour. 40 Figure 3.2 Converted F431A HCN2 data for activation and deactivation: Converted open probability during activation and deactivation for F431A HCN2 channel. Deactivation closing kinetic is alerted compared to wild type; it is much delayed with less voltage sensitivity. 41 3.1.3. Overlapped wild type and F431A HCN2 data In order to quickly compare and visualize the differences between wild type and F431A mutant HCN2 channels, overlapped open probability graphs are prepared in Figure 3.3. Figure 3.3 is simply the collection of Figures 3.1 and 3.2. 42 Figure 3.3 Overlapped curves from wild type and F431A HCN2 data: Overlapped open probability for activation and deactivation in both wild type and F431A HCN2 channel. Although activation is relatively less effected, it also shows major difference in opening kinetics but not as big as deactivation has. 43 3.2. The 4-state model 3.2.1. Analytical solution and its reproduced data In order to verify the derived and customized 4-state model, analytical solution was developed. As previously discussed in the introduction and methods, the analytical solution of the 4-state was calculated manually. The solution provides important initial conditions and exact expressions for the four states at steady-state. Table 3.1 summarizes the exact expression and the analytical solution for the state. The S and P terms are also explicitly expressed using the eight parameters obtained from the 4-state model. Those S and P terms are arbitrary choice of letter to simplify the long and complicated expression. With the given voltages, each term becomes simple number. The four states become the normalized open probability as the total sum of the four states must be 1. 44 Table 3.1 Analytical solution of the 4-state model States Expression Or Oa Cr Ca \u00CE\u00B1 \u00CE\u00B2 P1 \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 P2 \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B2 \u00CE\u00B2 P3 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 P4 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 P5 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 P6 \u00CE\u00B1 \u00CE\u00B1\u00CE\u00B2 \u00CE\u00B2 P7 \u00CE\u00B2 \u00CE\u00B2 S1 S2 45 In order to verify if the analytical solution is correctly calculated, Siegelbaum's original data (Chen et al., 2007) was employed for comparison. Figure 3.4 is the result of the reproduced simulation for both activation and deactivation. Although exact replication of experimental data was not possible on the simulations, actual overlapping of the two simulations on the same scale using printed paper also confirmed that the fits were reasonable and similar to those presented by Chen et al.. This further supports that the analytical solution is not only correct, but it also justifies and verifies the constructed 4-state model in MATLAB code with general solution in eigenvalues and eigenvectors 46 Figure 3.4 Reproduced simulations of Chen et al using analytical solution: Reproduced simulations using analytical solution of 4-state model to compare with Chen et al. The direct comparison between my and their simulations is not done here but it is confirmed that they are identical. This quickly validates the analytical solution of 4-state model. 47 3.2.2. Wild type HCN2 channel fitting Figure 3.5 is the fitting of wild type data according to the4-state model. The graph is produced using the set of parameters found by MATLAB. 48 Figure 3.5 Wild type HCN2 data fitting: Comparison between 4-state model and wild type data for activation and deactivation. Performance of 4-state model is not acceptable in both activation and deactivation. 49 3.3. The 6-state model The fitting results of the 4-state model were not fully satisfactory. Thus, the 6-state model was used for the fitting and compared to 4-state model below. 3.3.1. Wild type HCN2 data fitting Wild type fitting with the 6-state model is displayed in the Figure 3.6. It is clear that the 6-state model performs significantly better than the 4-state model, especially for the closing. 50 Figure 3.6 Fitting of the 6-state model for the wild type HCN2 data: Comparison between 6-state model and wild type data for activation and deactivation. Performance of 6- state model is much better in both activation and deactivation. 51 3.3.2. 4-state model vs. 6-state model It is very clear that the 6-state model is better than the 4-state model but the logical explanation for better performance of the 6-state model is required first. We knew that reasonable fitting of wild type data has to be done before fitting with the F431A mutant data. The proposed explanation is that extra states in the 6-state model allow more detailed and sophisticated changes over the course of time. Figure 3.7 shows the two possible pathways for the opening process. Figure 3.7 Two possible pathways in opening: Two opening pathways are highlighted. They are majorly used two pathways among many different ones due to the chosen parameter values. Parameters are carefully estimated in order to achieve the system of two opening pathways. The activation of wild type shows very dynamic changes even only with a 10mV difference in potential. The single opening pathway of the 4-state model is not enough to handle such diversity. In the case of deactivation, there is only one single closing pathway due to the determined magnitude of the parameters in the 6-state model as well. However, even one 52 extra state in that single pathway (O1) effectively handles the sigmoidicity in the closing. Figure 3.8 and 3.9 are continuous changes in each states from the 4- and 6-state model during deactivation at 0mV and +20mV respectively. 53 Figure 3.8 Changes in individual state at 0mV: Individual state analysis for wild type deactivation at 0mV for 4-state (top) and 6-state (bottom) model. From such an analysis, one is easily able to distinguish more used pathways and transiting states that majorly contributes the gating of channel. 54 Figure 3.9 Changes in individual state at +20mV: Individual state analysis for wild type deactivation at 20mV for 4-state (top) and 6-state (bottom) model. From such an analysis, one is easily able to distinguish more used pathways and transiting states that majorly contributes the gating of channel. 55 From both Figures, there is a significant accumulation of O1 state (black curve) in the 6-state model simulation. A similar curve is also produced by Or state (green curve) in the 4- state model simulation. This reveals that O1 of the 6-state model acts like Or of the 4-state model but at a much later stage of closing. This late involvement of the O1 makes it possible for the 6-state model has more flexibility to fit the experimental data. The individual states and their changes are further analysed in the discussion section. This brief explanation provides the logical validation of the 6-state model before fitting F431A mutant data. 3.3.3. F431A mutant data fitting As the 6-state model fits wild type data with more confidence, it is used to fit the F431A mutant channel data. Figure 3.10 and 3.11 are fittings for the activation and deactivation data respectively. 56 Figure 3.10 F431A activation fitting with the 6-state model: Comparison between 6- state model and F431A data for activation. Performance of 6-state model is still acceptable and reasonable for most of voltage values. However, there are big disagreement for -120mV and -110mV due to the poor performance of fitting and possible experimental error for those voltage values. The indigenous ion channels of oocytes are possible reasons. 57 Figure 3.11 Deactivation fitting with the 6-state model (-20mV): Comparison between 6-state model and F431A data for deactivation at -20mV. Performance of 6-state model is still acceptable as it well captured major features of closing kinetic of mutant channel. 58 Figure 3.12 Deactivation fitting with the 6-state model (0mV): Comparison between 6-state model and F431A data for deactivation at 0mV. Performance of 6-state model is still acceptable as it well captured major features of closing kinetic of mutant channel. 59 Figure 3.13 Deactivation fitting with the 6-state model (+20mV): Comparison between 6-state model and F431A data for deactivation at +20mV. Performance of 6-state model is still acceptable as it well captured major features of closing kinetic of mutant channel. 60 3.4. Goodness of the fit It is very difficult to tell the level of fitting by simply looking at the individual graph. As it is explained in the method section, the difference formula is used to compare and evaluate the goodness of the fit between the two models. Table 3.2, 3.3 and 3.4 are resultant calculations for the fittings of both models. The average differences are calculated using least square difference formulas as previously discussed in the method. Table 3.2 Wild type fitting evaluation for the 4-state model Activation Voltage (mV) Average differences -120 3.37% -110 2.60% -100 2.29% -90 3.83% -80 4.43% -70 4.37% -60 7.46% -50 5.35% Deactivation Voltage (mV) Average differences 0 6.01% +20 9.34% 61 Table 3.3 Wild type fitting evaluation for the 6-state model Activation Voltage (mV) Average differences -120 1.30% -110 1.69% -100 2.79% -90 2.80% -80 0.88% -70 1.94% -60 2.04% -50 2.79% Deactivation Voltage (mV) Average differences 0 0.62% +20 1.08% 62 Table 3.4 F431A fitting evaluation for the 6-state model Activation Voltage (mV) Average differences -120 4.59% -110 5.82% -100 1.47% -90 1.06% -80 1.34% -70 1.76% -60 1.46% -50 0.69% Deactivation Voltage (mV) Average differences -20 1.10% 0 1.20% +20 1.71% 63 3.5. Parameters Table 3.5, 3.6 and 3.7 are values of the parameters that were used for the 4-state and 6-state model fitting. For the 6-state model, different sets of parameters were used for wild type and the F431A mutant data. It is important to realize that the parameters are not just simple numbers. They all have physical meaning that is related to the gating of HCN channel. The unit for both voltage-independent and dependent parameter is as they indicate how fast is the transition among connected states. The exponent terms in the voltage- dependent parameter represents the voltage sensitivity and its direction of value change. terms increase while voltage become more negative and terms increase while voltage become more positive. The degree of such a change is determined by exponent terms. Table 3.5 Parameters for the 4-state model Voltage-independent parameters ( ) Voltage-dependent parameters ( ) 64 Table 3.6 Parameters for the 6-state model (wild type HCN2 channel) Parameters for wild type fitting Voltage-independent parameters( ) Voltage-dependent parameters ( ) \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 Same as \u00CE\u00B1 Same as \u00CE\u00B2 Same as \u00CE\u00B1' Same as \u00CE\u00B2' Table 3.7 Parameters for the 6-state model (F431A mutant HCN2 channel) Parameters for F431A mutant fitting Voltage-independent parameters ( ) Voltage-dependent parameters ( ) \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B1 \u00CE\u00B2 Same as \u00CE\u00B1 Same as \u00CE\u00B2 Same as \u00CE\u00B1' Same as \u00CE\u00B2' 65 4. Discussion 4.1. Motivation for the mathematical model Prior to analyzing the results, it is worthwhile to discuss the motivation behind mathematical modelling. For considerable amount of time, scientists and scholars have been investigating the world in order to fully understand it. There had been major and rapid discoveries, development, and growth in scientific technology and in the understanding of the universe from subatomic molecules to planetary objects. Each individual academic discipline provides fairly independent contribution towards the overall understanding of science for a long time. However, now it is common to study and research the same question with different academic disciplines. Hence, there are new fields such as biotechnology, biochemistry, biophysics, bioinformatics, biomaterials, and more. A classic yet relatively modern example is the discovery of the DNA\u00E2\u0080\u0099s double helix structure in 1952. It was the collaboration between James Watson and Francis Crick. James Watson is a biologist and Francis Crick is a physicist. These two very different academic backgrounds made a great synergy that led them to winning a Nobel Prize in 1962. 50 years later, the diverse collaboration between different disciplines to solve problems is expected and required as world is much more complicated than what was anticipated years ago. In order to understand and reveal the full complexity of the world, reasonable and careful estimations have to be made due to the limitations of our current knowledge and technology. Despite the many advanced experimental instruments and methods, there are still numerous cases in which researchers have to assume the underlying mechanisms from observations at a macroscopic scale. This is where mathematical models and simulations can 66 be very useful. Mathematical modelling is a powerful tool as long as reasonable assumptions are provided. Dynamic models are simplified representations of a particular real-space entity in the form of equations or computer code. They are intended to mimic essential features of the study system while leaving out nonessential features. The models are dynamic because they describe how system properties change over time. Mathematics, computation, and computer simulations are used to analyze models and to produce results that predict and describe the study system. The construction of successful models is constrained by what we can measure, either by estimating parameters that are part of the model formulation or by validating model predictions. Dynamic models are invaluable because they allow us to examine relationships that cannot be sorted out by purely experimental methods. A dynamic or mathematical model is built from equations representing the processes thought to account for the patterns observed in the data. The most powerful and useful property of a model is that it embodies a hypothesis about the study system and allows one to compare hypothesis with data without actual experiments. The model itself is built upon the observed patterns from the experimental data. Hopefully, if a model has sufficient characteristics of the study system to describe and mimic it, one can claim the original hypothesis is correct because the root of the model came from that hypothesis. 67 4.2. The 4-state model 4.2.1. Fitting of channel activation The fitting of the 4-state model for activation is yielded a range of average difference between 2.29% to 7.46%. The initial sigmoidal pattern of some curves show a reasonable level of fitting at some voltages, but not all. The biggest issue with the 4-state model with wild type activation is that it is really hard to fit every voltage value from -120mV to -50mV. That is, if certain set of parameter fits better at more hyperpolarizing voltages (-120mV,- 110mV-100mV), the fitting at less hyperpolarizing voltages (-90mV,-80mV,-70mV,-60mV,- 50mV) becomes poor. The exact opposite pattern was observed when the opposing voltages are favoured. The possible explanation for this issue is that the running time of simulation for each of the voltage value is quite different. Table 4.1 summarizes different running time. Table 4.1 Wild type experimental time Voltage (mV) Running time More hyperpolarizing group -120 10000 time units (2 seconds) -110 25000 time units (5 seconds) -100 25000 time units (5 seconds) Less hyperpolarizing group -90 50000 time units (10 seconds) -80 75000 time units (12 seconds) -70 75000 time units (12 seconds) -60 75000 time units (12 seconds) -50 75000 time units (12 seconds) 68 The length of running time of the shortest and the longest differs by 6 times. The more hyperpolarizing voltage group ran for a considerably shorter time frame compared to the other group. On top of that, the more hyperpolarizing voltage group had to be activated and open much faster than the other. Their kinetics is much faster even in a much shorter duration. These extreme behaviours from one group to another cannot be described nor handled by the simple 4-state model. 4.2.2. Fitting of the deactivation It is observed that the fitting of closing is worse than that of opening. The average differences of closing are 6.01% for 0mV and 9.34% for +20mV. The main difficulty with fitting the 4-state model is that the closing time of the simulations is either too fast or too slow compared to the experimental data. There is no median curve that can reasonably capture the behaviour of wild type deactivation. The fit of our data by the 4-state model was not as good as it was for the data from Chen at al. Why this is the case is not clear. One possibility is that Chen et al used various sets of parameters instead of one; from the table of parameters from the Chen et al, certain parameter values were given as a range. This suggests that actual fittings and simulations were done with various parameters in order to get the best fit. This would explain the difference between the original and our 4-statde model in the degree of the fitting, especially with the closing data. Another possibility is that our set of parameters may not be actually the best, although this would require a strict mathematical comparison. The model is sufficient as long as it describes the distinctive features of the target system. Lastly, some differences may be due to the fact that the experiments by Chen et al were carried out by patch clamping 69 in patches excised from the oocyte expressing HCN2, whereas our experiments were carried out using the two microelectrode voltage clamp. For example, it is possible that intracellular components in our experiments, such as cAMP and PIP2, might confer greater complexity such as sigmodicity in the recorded currents that are not well-captured by the 4 state model. Also, the range of activation is much less negative in our 2-micro-electrodes experiments. Our 6-state model HCN data in both wild type and F431A mutant channel fit better than the 4-state model. The extensive study and logical explanation to support the 6-state model were done to fully replace the 4-state model. 4.2.3. Analytical solution The analytical solution of the 4-state model serves well as the proper initial conditions. It also reproduces the original simulations to validate our own 4-state model. However, there are some limitations and considerations for the analytical solution despite its initial success. Initially, it was used to provide more insights on the parameters for the states since the analytical solution is the explicit expression in terms of parameters. Due to the complexity in the explicit expressions of states, they are not that useful even though a significant amount of time and effort were spent on manually solving the analytical solution. Although the analytical solution may still be valuable as definite initial conditions for the 6-state model, required initial conditions for activation and deactivation can be estimated. Since both protocols are from extreme ends of either hyperpolarizing or depolarizing voltage values, one can easily predict that most channels are initially open or closed. The written MATLAB code and constructed 4-state model were verified by reproducing original simulations and modifying them for the 6-state model is not completely new. Therefore, analytical solution of 70 the 6-state model is not needed nor does it have to be calculated as there are very limited advantages and information that can be drawn from it. 71 4.3. The 6-state model 4.3.1. Fitting of the wild type current data The fitting of the 6-state model is much better than that of the 4-state model. The biggest average difference is only 2.8% from between activation and deactivation. That value is actually close to the smallest average difference (2.29%) for the 4-state model. As previously stated, the biggest issue with the activation fitting is the diversity in opening kinetics for two different groups, namely the more (-120mV, -110mV, -100mV) and the less (-90mV,-80mV,-70mV,-60mV,-50mV) hyperpolarizing voltages groups. As it is shown in Figure 3.7, there are two possible pathways for the opening process in the 6-state model. These two different pathways would be utilized differently for the more and less hyperpolarizing voltage groups. The two voltage groups take different pathways: the first pathway as Cr\u00E2\u0086\u0092C1\u00E2\u0086\u0092O1\u00E2\u0086\u0092O2 and second pathway as Cr\u00E2\u0086\u0092C1\u00E2\u0086\u0092C2\u00E2\u0086\u0092O2. The only difference between the two pathways is the presence of O1 or C2. In order to validate which pathway is used, changes in individual state are plotted for the two voltage groups. Specific states show accumulation and increase in open probability changes (as states undergo a transient process for activation. Figure 4.1 and 4.2 are simulation result for the more and less hyperpolarizing voltage groups, respectively. The important curve to review is the solid blue line representing the population or open probability of C2 state. We know that the high accumulation of C2 is only possible if the opening process uses the second pathway. There is a strong accumulation of C2 for the more hyperpolarizing voltages group but insignificant accumulation of C2 for the less hyperpolarizing voltages group. This striking information reveals that the current elicited by more hyperpolarizing voltages, which open more quickly, takes the full advantage 72 of two pathways. They use both pathways at same time in order to achieve desirable level of open probability in shorter length of time. Meanwhile, the currents elicited by less hyperpolarizing voltages use mainly the first pathway and, thus, they have a much slower rate of opening and a lower level of open probability. They do not have to take the second pathway because their activation kinetic is not as rapid. Furthermore, more hyperpolarizing voltages are more likely inducing the movement of voltage sensors. Arriving to the C2 state requires movement of all four voltage sensors but less hyperpolarizing voltages group is not strong enough to cause the movement of all voltage sensors. This is exactly why the 4-state model, with only single pathway, is insufficient for fitting both groups. The similar individual state simulation is already provided as Figure 3.8 and 3.9. There is only single closing pathway for both the 4-state and 6-state model. However, there is very distinctive pattern is observed. The extra state of O1 plays a key role to adjust and enhances the overall fitting of the 6-state model. In order to properly describe the initial sigmoidal pattern, which is due to the slower and delayed kinetic at the beginning, O1 has to be increased to maintain the relatively high population of total open probability. After this initially slower phase, closing becomes accelerated so decrease in population of total open probability is needed. This can be managed by gradual decrease of O1 state until end of the simulation. This delicate control and management of O1 confirm that this state acts like a reservoir of the system. The similar task was given to the Or state for the 4-state model. However, without contribution of other extra states, the Or was not able to fit the wild type data properly. In the 4-state model, Or is only able to move to Cr state. The total open probability of the 4-state model cannot be finely managed as the channel has to be closed once it reaches Or state. That is why simulated closing of the 4-state model is either too slow 73 or too fast compared to the actual experimental data. Although, only the C2 and O1 are primarily emphasized from the 6-state model, it is important to realize that all of six states in the model are involved and contributed to the fitting. In other words, there are non-used states or meaningless states to be left. This concludes that the 6-state model is the simplest model to successfully describe the wild type of HCN channel. 74 Figure 4.1 More hyperpolarizing voltage group activation process for wild type HCN2 channel: Individual state analysis of 6-state model for the wild type activation. Solid blue curve for C2 state is only highly accumulated for -120mV and -110mV. There is still visible but small accumulation of C2 is observed for -100mV. 75 Figure 4.2 Less hyperpolarizing voltage group activation process for wild type HCN2 channel: Individual state analysis of 6-state model for the wild type activation. Solid blue curve for C2 state is not highly accumulated and not visible for -90mV, -80mV, and -70mV. This indicates at more positive potential, second pathway is not used as much. 76 4.3.2. The fitting of F431A mutant channel The fitting of F431A mutant channel by the 6-state model was acceptable as well, and better that the wild type channel as most average differences are in range of 1%. However, there are major discrepancies in the opening for 120mV (4.59%) and -110mV (5.82%). This may be due to the unstable reading and possible experimental error. Every open probability is normalized to -120mV where there is supposedly the greatest and maximal open probability. However, normalized open probability of -110mV is slightly great than 1 (around 1.05). The mutant channel also has more unstable holding potentials, especially at more hyperpolarizing voltages. This may indicate the unreliable experimental data for -120mV and -110mV. Same individual analysis was performed for F431A mutant data and it confirmed that the mechanism of gating for activation is similar between wild type and mutated channel. As shown in Figure 4.3 and 4.4, opening of F431A mutant channel still requires the two pathway system just like wild type channel does. It should be expected as there is no major difference in the activation curves. Since the rate of activation is similar, their mechanism should be similar as well. The deactivation process is similar even though mutated channel exhibits much delayed process. Significant and high involvement of O1 state further indicates the importance of this extra state along with the other five states. 77 Figure 4.3 More hyperpolarizing voltage group activation process for F431A: Individual state analysis of 6-state model for F431A mutant channel for the activation. Similar to wild type, two pathway opening system is still existed as clear accumulation of blue curve in -120mV,-110mV, and -100mV. 78 Figure 4.4 Less hyperpolarizing voltage group opening process for F431A: Individual state analysis of 6-state model for F431A mutant channel for the activation. Similar to wild type channel, there is no accumulation of C2 state for -90mV, -80mV, and - 70mV. This again confirms more positive (less hyperpolarizing) voltages do not use second pathway for the opening. 79 4.4. Parameters With the aid of MATALB and numerous simulations, two sets of the parameter were examined for the wild type and F431A mutant data. Originally, it is suspected that altering voltage independent parameter for closing (h' or g') should be enough to reasonably describe the F431A data, since mutant channel shows notable changes only during the deactivation. However, it turns out that many parameters have to be changed as shown in Figure 4.5. Red symbols indicate parameters with an increased value while blue symbols signify a decrease. Figure 4.5 Modified parameters: Summarized diagram to indicate modified parameters. Parameters in red are increased and ones in blue are decreased from wild type HCN2 channel. 4.4.1. Voltage dependent parameters Table 4.2 summarizes the percentage changes of the voltage dependent parameters compared to the wild type parameters. The positive percentage means that particular parameter has been increased by a given percentage from the value of the wild type channel and vice versa. 80 Table 4.2 Percentage changes for the voltage dependent parameters Voltage(mV) \u00CE\u00B1 \u00CE\u00B2 \u00CE\u00B2\u00E2\u0080\u0099 -120 96.657% 812.618% -72.299% -110 100.630% 704.175% -78.275% -100 104.683% 608.619% -82.962% -90 108.818% 524.416% -86.637% -80 113.036% 450.220% -89.520% -70 117.340% 384.839% -91.781% -60 121.730% 327.228% -93.554% -50 126.209% 276.462% -94.945% -40 130.779% 231.729% -96.035% -30 135.441% 192.311% -96.890% -20 140.197% 157.577% -97.561% -10 145.050% 126.970% -98.087% 0 150.000% 100.000% -98.500% 10 155.050% 76.235% -98.824% 20 160.203% 55.294% -99.077% 30 165.459% 36.841% -99.276% 40 170.822% 20.580% -99.433% It is quite surprising to see that voltage dependent parameters for the opening process are affected and their changes are very large as well. The reduction in \u00CE\u00B2' is expected for the delayed deactivation process. The large increases in \u00CE\u00B1 and \u00CE\u00B2 must be due to the fact 81 that F431A mutation induces changes in the overall opening process. Figure 3.3 shows the overlapping open probability for wild type and F431A mutant channels at every voltage value. Although the activation curves between wild type and mutant channels do not have statistically significant differences, the actual activation processes are quite different. This is why certain voltage dependent parameters in the activation process are modified significantly. It is also true that the modified \u00CE\u00B1 and \u00CE\u00B2 are not responsible for the delayed deactivation process as closing simulations were not significantly altered by them. 4.4.2. Voltage independent parameters Table 4.3 summarizes the percent changes of the voltage independent parameters. Table 4.3 Percentage changes in voltage independent parameters h g' h' h'' -25.000% -95.440% -99.544% 30.000% The near100% reduction in g' and h' are observed and they must be the primary and main contributor for the altered deactivation behaviour in the mutant channel. However, notable amount of reduction in the voltage dependent parameter (\u00CE\u00B2') is also found. One must confirm which parameter is truly responsible for the deactivation behaviour. A quick simulation had been performed with two different sets of parameter. Using our wild type channel data, both sets had identical wild type parameters except one set had 1000 times smaller b\u00E2\u0080\u0099 (voltage dependent parameter) and another set had 1000 times smaller g' and h' (voltage independent parameters). One set is expected to behave like the F431A mutant 82 channel, identifying the set that is responsible for delayed closing. If the model with modified g' and h' describes better than that of modified b\u00E2\u0080\u0099, then voltage independent parameters are more important than voltage dependent parameter for the deactivation of mutant channel. Figure 4.6 is the result of the simulation. Only the deactivation process was simulated as the purpose of this was to find more important parameters for the deactivation process. 83 Figure 4.6 Deactivation simulations of two models: Two group simulations for deactivation. The model with only reduced voltage-dependent parameters is only able to capture the delayed behaviour of mutant channel closing. But, mode with only reduced voltage-independent parameters is able to describe full features of mutant channel closing: delayed closing and much enhanced voltage independent behaviour. 84 The upper curves were produced by model with modified voltage dependent parameter and lower curves were produced by model with modified voltage independent parameters. It is very clear that modifying only voltage dependent parameter alone is not enough to accurately describe the true F431A behaviour. It does delay the deactivation time significantly compared to the wild type channel but there is no enhanced voltage independent behaviour. The deactivation of F431A channel is not only slow but it also has enhanced voltage independent behaviour. This finding eliminates the involvement of voltage dependent parameter as the main factor. The reproduced curves from the modified voltage independent parameters resemble closely to true F431A data. They show the delayed and strong voltage independent characteristics as some of the curves completely overlap with each others. This simulation clearly shows that voltage independent parameters are more important and are significant contributors to the deactivation behaviour of the F431A mutant channel. Moreover, reduced voltage dependent parameter for the closing is required to provide certain degree of voltage sensitivity to the mutant channel. As previously shown in Figure 3.2, more depolarizing voltage has faster deactivation kinetics like the wild type channel does but with less sensitivity. The combined reductions in both voltage dependent and independent parameters for closing are required to fully describe the mutant experimental data. However, the main responsible parameter for altered deactivation behaviour of the mutant channel, which was the initial motivation of the study, is indeed the voltage independent parameters. 85 4.5. Physiological significance The present study and analysis lead to a better understanding of the gating mechanism for HCN2 channel. The F431A mutation mainly alters the voltage-independent step by reducing its rate. Since the degree of reduction is so significant, the voltage- independent step becomes the rate determining step. The fact that the slowest step is the rate determining step is not new but the significance of the study is that this step specifically identifies the responsible location at the pore of HCN2 channel. This study is the first to suggest that the pore of the HCN2 channel controls voltage-independent opening. A voltage- independent transition is not only academically interesting, but it also has physiological meaning in terms of cardiac rhythmicity. HCN channel is known for its pacemaking activity and thus the discovery of a voltage-independent step can be used as a new target to control heart rate. With the two distinctive ways to activate and deactivate the channel, there will be more options and possibilities to regulate heart rate upon solid knowledge of this voltage- independent step. Two very well known drug blockers, Ivabradine and ZD7288, are proposed to bind inside of pore, where the responsible spot for voltage-independent step is located (BoSmith et al., 1993; Bucchi et al., 2002). They are also open-channel blockers. My simulations show that channel is always open regardless of voltage as low as 2% even in the wild type HCN2 channel. Furthermore, F431A mutation in HCN2 channel pore induces much delayed closing that drugs would have enough time to be utilized more efficiently and strongly regardless of applied voltage. 86 4.6. Limitations Although, the 6-state model was able to fit both the wild type and mutant channel data, there are possible limitations and issues presented. First, it is not certain that set of parameters we ultimately obtained is indeed the best. In fact, fitting for -120mV and -120mV of mutant activation data were not perfect. Generally, this would be one of most questionable issues from the constructing any mathematical model. Mathematical model can be a powerful tool for countless applications as long as valid assumptions are provided. It is like writing a program code to accomplish the same task. There will be multiple ways and approaches to complete same task in a different program coding. Similar to our mathematical model and its parameters, there will be multiple and different sets of parameters to describe same experimental data, although it does not seem likely that the overall level of fitting would be improved. Second, some of experimental data do not reach the true steady-state. For example, deactivation of F431A mutant channel data does not have a plateau region at the end of the experiments. This does not necessarily add the extra difficulty for the fitting but many of simulations tend to be in steady-state at the end. This raises some issues in degree of overall fitting as simulations are already in steady-state but actual experimental data are not. These disagreements may be overcome and improved as actual experimental data has true steady- state at the end. For future experiments, we can increase the running time to ensure every voltage value reaches the steady-state. On top of that, it will be better and interesting to see if every voltage value has the same longer running time because diverse running time for current activation protocol was first initial challenge for the data fitting. 87 4.7. Future applications and possible improvements The 6-state model can be used to describe different isoforms of HCN channels besides HCN2 channel. The 6-state model can be used to fit original HCN1 mutant data. It can also be used for HCN (homologs) channels from different species such as sea urchin or drosophila. Unfortunately, the 6-state model does not reflect the dual-allosteric of HCN channels because it does not have included effect of binding of cAMP. There are very active researches to fully understand true picture of cAMP binding including but not limited to its affinity, cooperativity, and required numbers of binding (Benndorf et al., 2012; Chow et al., 2012). Moreover, binding of different cyclic nucleotide is primarily studied in our lab. One of the biggest advantages of mathematical modelling is its ability to predict outcomes of experiments. Once the present 6-state model includes the binding of cAMP, it can be used to predict and compare the simulation with actual experimental data to strength and to support the findings. It will be also interesting to see if the binding of cAMP mainly influences only one type of steps or not between voltage-dependent and independent step. 88 5. Conclusion The delayed and enhanced voltage independent behavior in deactivation from F431A mutant channel can be explained by mainly modifying voltage independent parameters for deactivation process of the 6-state model. The 6-state model is the simplest model that reasonably and accurately describes both our experimental data of wild type and F431A mutant channel to support presence of voltage-independent step. 89 References Accili EA, Proenza C, Baruscotti M & DiFrancesco D. (2002). From funny current to HCN channels: 20 years of excitation. News Physiol Sci 17, 32-37. Altomare C, Bucchi A, Camatini E, Baruscotti M, Viscomi C, Moroni A & DiFrancesco D. (2001). Integrated allosteric model of voltage gating of HCN channels. J Gen Physiol 117, 519-532. Baruscotti M, Barbuti A & Bucchi A. (2010). The cardiac pacemaker current. J Mol Cell Cardiol 48, 55-64. Baruscotti M, Bucchi A & Difrancesco D. (2005). Physiology and pharmacology of the cardiac pacemaker (\"funny\") current. Pharmacol Ther 107, 59-79. Bell DC, Yao H, Saenger RC, Riley JH & Siegelbaum SA. (2004). Changes in local S4 environment provide a voltage-sensing mechanism for mammalian hyperpolarization-activated HCN channels. J Gen Physiol 123, 5-19. Benarroch EE. (2013). HCN channels: function and clinical implications. Neurology 80, 304-310. Benndorf K, Thon S & Schulz E. (2012). Unraveling subunit cooperativity in homotetrameric HCN2 channels. Biophys J 103, 1860-1869. Biel M, Wahl-Schott C, Michalakis S & Zong X. (2009). Hyperpolarization-activated cation channels: from genes to function. Physiol Rev 89, 847-885. BoSmith RE, Briggs I & Sturgess NC. (1993). Inhibitory actions of ZENECA ZD7288 on whole-cell hyperpolarization activated inward current (If) in guinea-pig dissociated sinoatrial node cells. Br J Pharmacol 110, 343-349. Bruening-Wright A & Larsson HP. (2007). Slow conformational changes of the voltage sensor during the mode shift in hyperpolarization-activated cyclic-nucleotide-gated channels. J Neurosci 27, 270-278. Bucchi A, Baruscotti M & DiFrancesco D. (2002). Current-dependent block of rabbit sino-atrial node I(f) channels by ivabradine. J Gen Physiol 120, 1-13. Chan YC, Wang K, Au KW, Lau CP, Tse HF & Li RA. (2009). Probing the bradycardic drug binding receptor of HCN-encoded pacemaker channels. Pflugers Arch 459, 25-38. 90 Chen J, Mitcheson JS, Tristani-Firouzi M, Lin M & Sanguinetti MC. (2001). The S4-S5 linker couples voltage sensing and activation of pacemaker channels. Proc Natl Acad Sci U S A 98, 11277-11282. Chen S, Wang J, Zhou L, George MS & Siegelbaum SA. (2007). Voltage sensor movement and cAMP binding allosterically regulate an inherently voltage-independent closed-open transition in HCN channels. J Gen Physiol 129, 175-188. Chow SS, Van Petegem F & Accili EA. (2012). Energetics of cyclic AMP binding to HCN channel C terminus reveal negative cooperativity. J Biol Chem 287, 600-606. Craven KB & Zagotta WN. (2004). Salt bridges and gating in the COOH-terminal region of HCN2 and CNGA1 channels. J Gen Physiol 124, 663-677. Craven KB & Zagotta WN. (2006). CNG and HCN channels: two peas, one pod. Annu Rev Physiol 68, 375-401. DiFrancesco D. (1984). Characterization of the pace-maker current kinetics in calf Purkinje fibres. J Physiol 348, 341-367. DiFrancesco D. (1993). Pacemaker mechanisms in cardiac tissue. Annu Rev Physiol 55, 455-472. DiFrancesco D. (1999). Dual allosteric modulation of pacemaker (f) channels by cAMP and voltage in rabbit SA node. J Physiol 515 ( Pt 2), 367-376. DiFrancesco D & Tortora P. (1991). Direct activation of cardiac pacemaker channels by intracellular cyclic AMP. Nature 351, 145-147. Gauss R, Seifert R & Kaupp UB. (1998). Molecular identification of a hyperpolarization-activated channel in sea urchin sperm. Nature 393, 583-587. Ishii TM, Takano M, Xie LH, Noma A & Ohmori H. (1999). Molecular characterization of the hyperpolarization-activated cation channel in rabbit heart sinoatrial node. J Biol Chem 274, 12835-12839. Lakatta EG & DiFrancesco D. (2009). What keeps us ticking: a funny current, a calcium clock, or both? J Mol Cell Cardiol 47, 157-170. 91 Ludwig A, Zong X, Jeglitsch M, Hofmann F & Biel M. (1998). A family of hyperpolarization- activated mammalian cation channels. Nature 393, 587-591. Macri V & Accili EA. (2004). Structural elements of instantaneous and slow gating in hyperpolarization-activated cyclic nucleotide-gated channels. J Biol Chem 279, 16832-16846. Mannikko R, Elinder F & Larsson HP. (2002). Voltage-sensing mechanism is conserved among ion channels gated by opposite voltages. Nature 419, 837-841. Mannikko R, Pandey S, Larsson HP & Elinder F. (2005). Hysteresis in the voltage dependence of HCN channels: conversion between two modes affects pacemaker properties. J Gen Physiol 125, 305-326. Postea O & Biel M. (2011). Exploring HCN channels as novel drug targets. Nat Rev Drug Discov 10, 903-914. Proenza C, Angoli D, Agranovich E, Macri V & Accili EA. (2002). Pacemaker channels produce an instantaneous current. J Biol Chem 277, 5101-5109. Proenza C & Yellen G. (2006). Distinct populations of HCN pacemaker channels produce voltage- dependent and voltage-independent currents. J Gen Physiol 127, 183-190. Ryu S & Yellen G. (2012). Charge movement in gating-locked HCN channels reveals weak coupling of voltage sensors and gate. J Gen Physiol 140, 469-479. Santoro B, Liu DT, Yao H, Bartsch D, Kandel ER, Siegelbaum SA & Tibbs GR. (1998). Identification of a gene encoding a hyperpolarization-activated pacemaker channel of brain. Cell 93, 717-729. Sherwood L & Cengage Learning (Firm). (2010). Human physiology : from cells to systems. Brooks/Cole, Cengage Learning, Australia ; United States. Vemana S, Pandey S & Larsson HP. (2004). S4 movement in a mammalian HCN channel. J Gen Physiol 123, 21-32. Wainger BJ, DeGennaro M, Santoro B, Siegelbaum SA & Tibbs GR. (2001). Molecular mechanism of cAMP modulation of HCN pacemaker channels. Nature 411, 805-810. 92 Appendix Appendix A : MATLAB code 4-state model simulation with analytical solution v=-105; g1=0.0024; h1=0.00038; g2=0.0000042; h2=0.021; a0=0.0000032; sa=9.1; b0=415.8; sb=49; c0=0.000011; sc=9.9; d0=0.045; sd=33.9; %Define the voltage-dependent parameters, which are also come from %Sigelbaum's paper. a1=a0*exp(-v/sa); b1=b0*exp(v/sb); a2=c0*exp(-v/sc); b2=d0*exp(v/sd); %B=[-a1 b1 0 0;a1 -(b1+g1) h1 0;0 g1 -(b2+h1) a2;0 0 b2 -a2]; %C=[-(a1) b1 0 0;a1 -(b1+g1) h1 0;0 g1 -(h1) a2;0 0 0 -(a2)]; A=[-(a1+g2) b1 0 h2;a1 -(b1+g1) h1 0;0 g1 -(b2+h1) a2;g2 0 b2 -(a2+h2)]; %A=[-(a1) b1 0 0;a1 -(b1+g1) 0 0;0 g1 -(0) a2;0 0 0 -(a2+0)]; [V,D]=eig(A); Eigenvalues=[D(1,1);D(2,2);D(3,3);D(4,4)]; %Getting ICs v=0; g1=0.0024; h1=0.00038; g2=0.0000042; h2=0.021; a0=0.0000032; sa=9.1; b0=415.8; sb=49; c0=0.000011; sc=9.9; d0=0.045; sd=33.9; 93 %Define the voltage-dependent parameters, which are also come from %Sigelbaum's paper. a1=a0*exp(-v/sa); b1=b0*exp(v/sb); a2=c0*exp(-v/sc); b2=d0*exp(v/sd); %Define expression for each of 4 states P4=(b1+g1)/(a1+b1+g1); P5=(b1+h1+g1)/(a1+b1+g1); P1=((a2+h2)+g2*(b1+g1)/(a1+b1+g1)); P2=((b2-g2*(b1+h1+g1)/(a1+b1+g1))^-1); P3=(g2*(b1+g1)/(a1+b1+g1)); P6=((a1+g2)/h2)-((b1*a1)/(h2*(b1+g1))); P7=b1*h1/(h2*(b1+g1)); S1=P4*P6+P2*P3*P5*P6-P2*P3*P7; S2=1-1*(-P4*P6-P1*P2*P5*P6-P1*P2*P7); %positive Or=S1/S2; Oa=P1*P2*Or-P3*P2; Cr=P4-P5*Oa-P4*Or; Ca=(a1*Cr+h1*Oa)/(b1+g1); tt=Or+Oa+Cr+Ca; ot=Or+Oa; ct=Cr+Ca; %Give initial condition at certain voltage value ICs=[Cr;Ca;Oa;Or]; Cs=V\ICs; x=0:1:8000; Cr_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(1,1)+Cs(2,1)*exp(Eigenvalues(2,1)*x )*V(1,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(1,3)+Cs(4,1)*exp(Eigenvalues(4, 1)*x)*V(1,4); Ca_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(2,1)+Cs(2,1)*exp(Eigenvalues(2,1)*x )*V(2,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(2,3)+Cs(4,1)*exp(Eigenvalues(4, 1)*x)*V(2,4); Oa_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(3,1)+Cs(2,1)*exp(Eigenvalues(2,1)*x )*V(3,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(3,3)+Cs(4,1)*exp(Eigenvalues(4, 1)*x)*V(3,4); Or_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(4,1)+Cs(2,1)*exp(Eigenvalues(2,1)*x )*V(4,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(4,3)+Cs(4,1)*exp(Eigenvalues(4, 1)*x)*V(4,4); C_t=Cr_g+Ca_g; O_t=Oa_g+Or_g; Max=0.8381; C_N=C_t./Max; O_N=O_t./Max; 94 %New=O_N.^2; %Figure(3) plot(x,O_N,'k') Codes for G.A. Initial response (customized MATLAB function) function y = initial_response20(x,t) x = x/10000; % Scale Down v = -60; %ac=4*10^-7; acc=0.075; apc=2.5*10^-6; apcc=0.045; bc=0.002; bcc=0.00115; bpc=0.0003; bpcc=0.0336; a = x(1)*exp(-x(2)*v); ap = apc*exp(-apcc*v); l = a; lp = ap; b = bc*exp(bcc*v); bp = bpc*exp(bpcc*v); d = b; dp = bp; g = 5*10^-4; gp = 0.57*10^-6; gpp =20*10^-4; h = 1.5*10^-4; hp = 0.57*10^-4; hpp = 0.65*10^-4; A = [ -(a+gp) b 0 0 0 hp; a -(b+l+g) d 0 h 0; 0 l -(d+gpp) hpp 0 0; 0 0 gpp -(dp+hpp) lp 0; 0 g 0 dp -(bp+lp+h) ap; gp 0 0 0 bp -(ap+hp) ]; B = [0 0 0 0 0 0]'; C = [1 1 1 0 0 0; 0 0 0 1 1 1]; D = 0; ss_system = ss(A,B,C,D); %y = initial(ss_system,[0 0 0 1 0 0],t); y = initial(ss_system,[1 0 0 0 0 0],t); dol_myFit (customized MATLAB function) function z = dol_myFit20(x,expData,t) y = initial_response20(x,t); penalty = sum((sum(y,2)-1).^2)*10000; z = sum((expData - y(:,2)).^2) + penalty; dol_GA (code for setting and running G.A.) %% Ctreating Dummy Data %t = (0:0.1:20)'; 95 %expData = 1 - exp(-t); t = 1:length(expData); %% GA Part options = gaoptimset; options.Generations = 500; %options.InitialPopulation = ones(1,8)*0.0001; nvars = 2; % Number of input Variables fitness = @(x) dol_myFit20(x,expData,t); lb = zeros(1,2); %lower bound for parameters ub = lb+inf; %upper bound for parameters [x,fval,exitflag,output,population,score] = ... ga(fitness,nvars,[],[],[],[],lb,ub,[],options); %x = x/10000; % Scale down %% Post Analysis y = initial_response20(x,t); plot(t,[expData y]); %axis([0 20 0 1.3]); xlabel('time(sec)'); legend('Experiment(Oa2+Oa1+Or)','Ca2+Ca1+Cr','Oa2+Oa1+Or'); 6-state simulation for individual analysis v = -120; x(1)=1*10^-6; x(2)=0.063; x(3)=2.5*10^-6; x(4)=0.045; x(5)=0.002; x(6)=0.00115; x(7)=0.00003; x(8)=0.0336; a = x(1)*exp(-x(2)*v); ap = x(3)*exp(-x(4)*v); l = a; lp = ap; b = x(5)*exp(x(6)*v); bp = x(7)*exp(x(8)*v); d = b; dp = bp; g = 5*10^-4; gp = 0.57*10^-6; gpp =20*10^-4; h = 1.5*10^-4; hp = 0.57*10^-4; hpp = 0.65*10^-4; A = [ -(a+gp) b 0 0 0 hp; a -(b+l+g) d 0 h 0; 0 l -(d+gpp) hpp 0 0; 0 0 gpp -(dp+hpp) lp 0; 0 g 0 dp -(bp+lp+h) ap; gp 0 0 0 bp -(ap+hp) ]; %B=[-a1 b1 0 0;a1 -(b1+g1) h1 0;0 g1 -(b2+h1) a2;0 0 b2 -a2]; %C=[-(a1) b1 0 0;a1 -(b1+g1) h1 0;0 g1 -(h1) a2;0 0 0 -(a2)]; %A=[-(a1+g2) b1 0 h2;a1 -(b1+g1) h1 0;0 g1 -(b2+h1) a2;g2 0 b2 -(a2+h2)]; %A=[-(a1) b1 0 0;a1 -(b1+g1) 0 0;0 g1 -(0) a2;0 0 0 -(a2+0)]; [V,D]=eig(A); Eigenvalues=[D(1,1);D(2,2);D(3,3);D(4,4);D(5,5);D(6,6)]; 96 %Getting ICs %Give initial condition at certain voltage value ICs=[1;0;0;0;0;0]; %ICs=[0;0;0;1;0;0]; %ICs=[Cr;Ca1;Ca2;Oa2;Oa1;Or]; Cs=V\ICs; x=0:75000; Cr_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(1,1)+Cs(2,1)*exp(Eigenvalues(2,1)*x )*V(1,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(1,3)+Cs(4,1)*exp(Eigenvalues(4, 1)*x)*V(1,4)+Cs(5,1)*exp(Eigenvalues(5,1)*x)*V(1,5)+Cs(6,1)*exp(Eigenvalue s(6,1)*x)*V(1,6); Ca1_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(2,1)+Cs(2,1)*exp(Eigenvalues(2,1)* x)*V(2,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(2,3)+Cs(4,1)*exp(Eigenvalues(4 ,1)*x)*V(2,4)+Cs(5,1)*exp(Eigenvalues(5,1)*x)*V(2,5)+Cs(6,1)*exp(Eigenvalu es(6,1)*x)*V(2,6); Ca2_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(3,1)+Cs(2,1)*exp(Eigenvalues(2,1)* x)*V(3,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(3,3)+Cs(4,1)*exp(Eigenvalues(4 ,1)*x)*V(3,4)+Cs(5,1)*exp(Eigenvalues(5,1)*x)*V(3,5)+Cs(6,1)*exp(Eigenvalu es(6,1)*x)*V(3,6); Oa2_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(4,1)+Cs(2,1)*exp(Eigenvalues(2,1)* x)*V(4,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(4,3)+Cs(4,1)*exp(Eigenvalues(4 ,1)*x)*V(4,4)+Cs(5,1)*exp(Eigenvalues(5,1)*x)*V(4,5)+Cs(6,1)*exp(Eigenvalu es(6,1)*x)*V(4,6); Oa1_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(5,1)+Cs(2,1)*exp(Eigenvalues(2,1)* x)*V(5,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(5,3)+Cs(4,1)*exp(Eigenvalues(4 ,1)*x)*V(5,4)+Cs(5,1)*exp(Eigenvalues(5,1)*x)*V(5,5)+Cs(6,1)*exp(Eigenvalu es(6,1)*x)*V(5,6); Or_g=Cs(1,1)*exp(Eigenvalues(1,1)*x)*V(6,1)+Cs(2,1)*exp(Eigenvalues(2,1)*x )*V(6,2)+Cs(3,1)*exp(Eigenvalues(3,1)*x)*V(6,3)+Cs(4,1)*exp(Eigenvalues(4, 1)*x)*V(6,4)+Cs(5,1)*exp(Eigenvalues(5,1)*x)*V(6,5)+Cs(6,1)*exp(Eigenvalue s(6,1)*x)*V(6,6); C_t=real(Cr_g+Ca1_g+Ca2_g); O_t=real(Oa2_g+Oa1_g+Or_g); %subplot(6,1,1) %plot(x,Cr_g) %subplot(6,1,2) %plot(x,Ca1_g) %subplot(6,1,3) %plot(x,Ca2_g) %subplot(6,1,4) %plot(x,Oa2_g) %subplot(6,1,5) %plot(x,Oa1_g) %subplot(6,1,6) %plot(x,Or_g) Figure(2) 97 if v==-120 plot(expData_a) end if v==-110 plot(expData2a) end if v==-100 plot(expData3a) end if v==-90 plot(expData4a) end if v==-80 plot(expData5a) end if v==-70 plot(expData6a) end if v==-60 plot(expData7a) end if v==-50 plot(expData8a) end if v==-20 plot(A1) end if v==0 plot(B1) end if v==10 plot(B1) end if v==20 plot(C1) end if v==40 plot(D1) end 98 hold on %Figure(2) plot(O_t,'r') "@en . "Thesis/Dissertation"@en . "2013-11"@en . "10.14288/1.0073987"@en . "eng"@en . "Cell and Developmental Biology"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "Attribution-NonCommercial-NoDerivatives 4.0 International"@en . "http://creativecommons.org/licenses/by-nc-nd/4.0/"@en . "Graduate"@en . "A kinetic model for the HCN2 ion channel"@en . "Text"@en . "http://hdl.handle.net/2429/44696"@en .