On Numerical Models of Calcium-Induced Calcium Release by Andrew T. L. Nixon B.S. Mathematics, Rensselaer Polytechnic Institute; Troy, NY, 2008 B.S. Physics, Rensselaer Polytechnic Institute; Troy, NY, 2008 M.S. Applied Mathematics, Rensselaer Polytechnic Institute; Troy, NY, 2008 A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The Division of Applied Mathematics at Brown University PROVIDENCE, RHODE ISLAND May 2015 c Copyright 2015 by Andrew T. L. Nixon This dissertation by Andrew T. L. Nixon is accepted in its present form by The Division of Applied Mathematics as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Bj¨orn Sandstede, Ph.D., Advisor Recommended to the Graduate Council Date Anastasios Matzavinos, Ph.D., Reader Date Martin Maxey, Ph.D., Reader Approved by the Graduate Council Date Peter M. Weber, Dean of the Graduate School iii Vitae Professional Preparation Rensselaer Polytechnic Institute Mathematics B.S., summa cum laude, 2008 Rensselaer Polytechnic Institute Physics B.S., summa cum laude, 2008 Rensselaer Polytechnic Institute Applied Mathematics M.S., 2008 Brown University Applied Mathematics Ph.D., expected 2015 Appointments Teaching Assistant, Rensselaer Polytechnic Institute Fall 2007 - Spring 2008 Research Assistant, Brown University Fall 2010 - present Teaching Assistant, Brown University Fall 2011 - Spring 2012, Fall 2013 - present Publications M. Crosskey, A. Nixon, L. Schick, G. Kovaˇciˇc. Invisibility cloaking via non-smooth transformation optics and ray tracing. Physics Letters A, 375, 1903-1911 (2011). iv Presentations Invisibility cloaks for asymmetric objects. SIAM Annual Meeting, Denver, 2009. Teaching Experience and Training MATH 1010-1020: Calculus I + II, Teaching Assistant, Rensselaer Polytechnic Institute APMA 0330-0360: Methods of Applied Math I + II, Teaching Assistant, Brown University MATH 0090: Introductory Calculus, Part I; Teaching Assistant, Brown University v Preface and Acknowledgments To Emma, for believing in me when I did not believe in myself. To my advisor, Bj¨orn Sandstede, for continuing to guide me and work with me through very difficult times. To my family and friends and everyone else who offered support and encourage- ment along this journey. vi Contents Vitae iv Preface and Acknowledgments vi 1 Introduction 1 1.1 The biological role of calcium . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Calcium ion channels . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Channel geometry and dynamics . . . . . . . . . . . . . . . . . . . . 7 1.4 At the channel: channel stochastics . . . . . . . . . . . . . . . . . . . 10 1.5 Between channels: reaction-diffusion . . . . . . . . . . . . . . . . . . 13 1.6 Coupling the regimes: the Gillespie algorithm . . . . . . . . . . . . . 17 1.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2 Approximations 22 2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2 Instantaneous equilibration . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.2 Fast timescale of channel influx . . . . . . . . . . . . . . . . . 27 2.2.3 Limited spatial extent of channel influx . . . . . . . . . . . . . 31 2.2.4 Eliminating time evolution from the diffusion computation . . 33 2.2.5 Handling the geometry . . . . . . . . . . . . . . . . . . . . . . 36 2.2.6 The Gillespie-R¨ udiger algorithm . . . . . . . . . . . . . . . . . 38 2.3 Constant buffers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.3.2 Justification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.4 Channel flux distribution . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.4.2 Localized flux . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.4.3 Uniform flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 2.5 Puffs as the fundamental event . . . . . . . . . . . . . . . . . . . . . 51 2.5.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.5.2 Computation of the cluster transition rates . . . . . . . . . . . 53 2.5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 vii 3 Channel Reaction Rates 58 3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2 Illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.3 Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4 Conclusion 74 4.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.2 Full model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2.1 Single channel . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2.2 Single cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2.3 The cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.3 Further work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 A The Gillespie-R¨ udiger Algorithm 82 A.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 A.2 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 A.3 The probability distribution . . . . . . . . . . . . . . . . . . . . . . . 88 A.4 Sampling the variables: Gillespie’s algorithm . . . . . . . . . . . . . . 90 A.5 Time dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 A.6 Solving the integral: R¨ udiger’s algorithm . . . . . . . . . . . . . . . . 97 B Analytic Steady-State Solution due to a Single Open Channel 99 B.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 B.2 Setup and reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 B.3 Flux condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 viii List of Tables 1.1 Channel transition rates from (De Young and Keizer, 1992). Refer to Figure 1.6 for a schematic illustrating the behavior of channel transi- tions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1 Comparison of channel rate constants from various models. Units are µM−1 s−1 for binding rates (ai ) and s−1 for dissociation rates (bi ). Repeated and elaborated upon in Chapter 3 (Table 3.2). . . . . . . . 30 3.1 Channel transition rates from Skupin et al. (2010). Refer to Figure 3.1 for a schematic illustrating the behavior of channel transitions. . . 60 3.2 Comparison of channel rate constants from various models: (a) Skupin et al. (2010), (b) De Young and Keizer (1992), (c) R¨ udiger et al. (2007), (d) Thurley and Falcke (2011), (e) those computed in this chapter. As before (Tables 1.1, 3.1), the ai are binding rates and have units µM−1 s−1 and the bi are dissociation rates with units s−1 . Note: many articles list values of ai and di where di ≡ abii . For these articles, we have instead computed and listed the values of bi for easier comparison. We have included the corresponding nomenclature used in this chapter for convenience. . . . . . . . . . . . . . . . . . . . . . 62 4.1 Summary of channel transition rates derived in Chapter 3. . . . . . . 77 ix List of Figures 1.1 Schematic of the role of calcium-induced calcium release in apoptosis. Compare with Figure 1.2. EGTA and BAPTA are examples of buffers (see Section 1.5). Figure recreated from Sanson et al. (2008). . . . . . 2 1.2 Schematic of the role of calcium-induced calcium release in cardiac excitation-contraction coupling. See text for details. Figure recreated from Viola and Hool (2010). . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 A cartoon depiction of different calcium channel events. Figure recre- ated from (Skupin et al., 2010). Note that blips and puffs are typically measured locally, whereas spikes are averaged over the volume of the cell, creating an apparent discrepancy in the concentration scales. . . 5 1.4 3D model of the protein structure of a voltage-gated potassium chan- nel. Note the four-fold symmetry. Figure recreated from Jiang et al. (2003). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.5 The fraction of time a single channel is open as a function of local calcium concentration. Includes both experimental data from (Mak et al., 1998) and simulation results from (R¨ udiger et al., 2007). Figure recreated from (R¨ udiger et al., 2007). . . . . . . . . . . . . . . . . . . 8 1.6 The original De Young-Keizer model. The quantities aj are binding rates and the bj are dissociation rates (Table 1.1). Figure recreated from (De Young and Keizer, 1992). . . . . . . . . . . . . . . . . . . . 10 2.1 Time evolution of the calcium concentration at an open channel with no approximations. Biological parameters used are from (R¨ udiger et al., 2007). Modeled using a finite-difference scheme with the full reaction-diffusion equations for c and bi with three species of buffer. Channel flux distribution uniform over the channel volume. . . . . . . 27 2.2 Time evolution of the spatial calcium concentration profile. Recreated from (Smith et al., 1996). Of particular note is the 10 mM exogenous buffer concentration. . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3 Time for the concentration at a point to reach 99.9% of the equilibrium concentration at that point. All model parameters are the same as in Figure 2.1 except the buffer species concentrations are treated as constant - both uniform throughout the cell and independent of time. The reason the time falls to 0 at around 1.3 µm is because c0 is already within 99.9% of the equilibrium value. Beyond this point, the effect of the channel is negligible. . . . . . . . . . . . . . . . . . . . . . . . . 29 2.4 Spatial distribution of calcium due to a cluster with five open channels 2.5 ms after opening. . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 x 2.5 Comparison of the spatial profiles due to the influx of a single chan- nel with (a) constant and (b) dynamically evolved buffers. Elapsed time is 2.5 ms. Biological parameters from R¨ udiger et al. (2007), with uniform flux over channel volume. Computed with 2nd-order radi- ally symmetric finite difference scheme and 4th-order Runge-Kutta method. Dirichlet conditions at cell boundary. See Figure 2.6 for a plot of the absolute error. . . . . . . . . . . . . . . . . . . . . . . . . 41 2.6 The absolute error corresponding to Figure 2.5. The error is small compared to c0 = .05 µM far from the channel. . . . . . . . . . . . . . 42 2.7 The time evolution of the buffer species concentrations at an open channel. Total elapsed time 2.5 ms. . . . . . . . . . . . . . . . . . . . 44 2.8 Steady state-concentration due to influx of a single open channel with the flux (a) spread uniformly over the channel volume, (b) localized at the channel center. For the localized source, the concentration profile is truncated at the channel boundary. . . . . . . . . . . . . . . . . . . 48 2.9 Steady state profile due to an open channel with uniform flux using a spatial granularity of (a) Rch , (b) R2ch , (c) R5ch , (d) R10ch . . . . . . . . . 49 2.10 A comparison of the resulting dynamics modeling a cluster as (a) con- sisting of 25 channels and using the traditional De Young-Keizer model, (b) a single stochastic object with transition rates as com- puted in Section 2.5.2. Note: the steep diagonal lines in part (b) result from automatic linear interpolation between two distant data points; the concentration is actually c0 on these intervals. . . . . . . . 55 2.11 Steady-state spatial profile due to a 5-channel-strength source (typical puff size) with the flux spread uniformly over (a) the volume of the cluster and (b) the volume of a single channel. The radius of a cluster is 30 nm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.12 A spike. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.1 The original De Young-Keizer model. The quantities aj are bind- ing rates and the bj are dissociation rates. Figure recreated from (De Young and Keizer, 1992). (Repeated from Chapter 1.) . . . . . . 61 3.2 Dynamics of a single cluster with transition rates for tetrameric chan- nels derived from experimental data. . . . . . . . . . . . . . . . . . . 72 3.3 Dynamics of a single cluster with parameters from (Swillens et al., 1999). Repetition of Figure 2.10a. . . . . . . . . . . . . . . . . . . . . 73 4.1 The fraction of time a single channel is open as a function of local calcium concentration. Repetition of Figure 1.5. . . . . . . . . . . . . 78 4.2 A comparison of the spatial equilibrium profile for (a) a single chan- nel with point source or uniform flux, (b) a cluster of channels with uniform flux over the volume of a channel vs. volume of a cluster. Repetition of Figures 2.8 and 2.11. . . . . . . . . . . . . . . . . . . . 79 xi Chapter ONE Introduction 2 Figure 1.1: Schematic of the role of calcium-induced calcium release in apoptosis. Compare with Figure 1.2. EGTA and BAPTA are examples of buffers (see Sec- tion 1.5). Figure recreated from Sanson et al. (2008). 1.1 The biological role of calcium Calcium is an important second messenger responsible for a variety of cellular and tissue processes including gene expression, cell division and differentiation, apop- tosis (Fig. 1.1), and muscle contraction (Fig. 1.2). In this work, we will study in detail models of cellular calcium dynamics. To provide context, however, we will first briefly outline one such process in which calcium is strongly implicated: car- diac excitation-contraction coupling (cardiac muscle contraction). This background will serve to illustrate the biological importance of calcium and help to motivate its study. In excitation-contraction coupling (Fig. 1.2) an action potential (originating from a motor neuron in skeletal muscle or from the sinoatrial node in the heart) travels through structures in the tissue called T-tubules. The action potential en- counters voltage-sensitive calcium channels called L-type calcium channels (LTCC), 3 Figure 1.2: Schematic of the role of calcium-induced calcium release in cardiac excitation-contraction coupling. See text for details. Figure recreated from Viola and Hool (2010). causing them to undergo a conformational change and open, allowing calcium ions to enter the cell cytoplasm. This calcium can then bind to IP3 receptors (IP3 R) or Ryanodine receptors (RyR) on channels at the surface of the sarcoplasmic reticu- lum (SR)1 . The binding of calcium to these channels induces them to open, releasing stores of yet more calcium ions from the interior of the SR into the cytoplasm, significantly amplifying the calcium signal. This process of amplification is called calcium-induced calcium release (CICR) and is the central topic of this thesis. In muscle tissue, this massive influx of calcium induces contraction of the cell (and thereby, the tissue) via a mechanism called the cross-bridge cycle. Long struc- tures called myofibrils stretch across a muscle cell. These structures are bundles of two different kinds of filaments: thinner filaments comprised of the protein actin, and thicker ones comprised of the protein myosin. Ca2+ ions bind to troponin pro- teins on an actin filament, causing them to undergo a conformational change and 1 The sarcoplasmic reticulum is a particular kind of endoplasmic reticulum (ER) and the two terms are used interchangeably in the context of calcium-induced calcium release. 4 expose binding sites on the actin filament. Protrusions on an adjacent myosin fil- ament (called myosin heads) can then attach to these binding sites. The myosin heads attach, pull on the actin filament, detach, and return to their original posi- tion. They do this repeatedly in a sort of ratcheting action, causing the actin and myosin filaments to slide against one another. When this occurs in all the myofibrils in the cell, it causes the cell to contract, and when all of the cells in a tissue undergo this process in concert, the tissue contracts. This illustrates the role calcium plays in translating an electrical impulse into the mechanical motion of a muscle tissue. In cardiac tissue, calcium-induced calcium release plays a critical role by amplifying the calcium signal. A fully quantitative understanding of the mechanisms driving calcium-induced calcium release can therefore provide a more complete understanding of excitation- contraction coupling, and in turn a more complete understanding of cardiac contrac- tion, creating new avenues for medical research and applications. In this work, we will therefore study individual calcium channels, study the mech- anisms that drive their behavior, and examine and improve upon quantitative nu- merical models of these channels and their behavior. This chapter will elaborate on the physical process of calcium-induced calcium release, including the importance of individual channels as well as the relevant geome- tries and dynamics, and will outline the details of the general approach to modeling CICR. Chapter 2 will discuss several approximations and simplifications we have investigated. Chapter 3 will examine the rate constants for channel binding and dissociation reactions, and provide a novel construction of a set of values for these rates determined from experimental data. Finally, Chapter 4 will summarize the results of Chapters 2 and 3, outline the full model as we have constructed it, and discuss further questions not answered by this work. 5 1.2 Calcium ion channels In the previous section, we briefly il- lustrated the biological role of calcium- induced calcium release, and how it is driven by the behavior of individual ion channels: in various kinds of cells, clusters of channels sensitive to inosi- tol 1,4,5-trisphosphate (IP3 ) or to Ryan- odine are responsible for releasing Ca2+ into the cytoplasm from reservoirs in the endoplasmic reticulum (ER). These channels are themselves gated by Ca2+ , creating a cooperative phenomenon; the random opening of one channel can en- courage other nearby channels to open Figure 1.3: A cartoon depiction of differ- ent calcium channel events. Figure recre- simultaneously. ated from (Skupin et al., 2010). Note that blips and puffs are typically measured lo- This cooperativity is significant be- cally, whereas spikes are averaged over the cause it creates a strongly positive feed- volume of the cell, creating an apparent discrepancy in the concentration scales. back process. Not only does this allow for the generation of complex signals, but it means that the termination of CICR must happen at the channel level, further motivating the study of the behavior of individual ion channels. Because of this cooperative behavior, we observe three types of calcium events (Fig. 1.3): blips - in which a single channel opens and closes repeatedly in rapid succession, puffs - which involve the concerted opening of several channels in one cluster, and spikes (sometimes “sparks”) - where several clusters in the cell are ac- 6 tive simultaneously. The complexity of these dynamics together with the wide variety of processes that they drive means that CICR is an area of study of significant interest, and many models have been created attempting to explain experimental data. This the- sis will discuss recent attempts to model individual calcium channels and clusters thereof. We will discuss in detail several existing approaches to modeling CICR and the assumptions used in these models. We will explore and develop a variety of new approximations and other assumptions and discuss the implications they have for the accuracy and efficiency of the model. 7 Figure 1.4: 3D model of the protein structure of a voltage-gated potassium channel. Note the four-fold symmetry. Figure recreated from Jiang et al. (2003). 1.3 Channel geometry and dynamics Imaging of ion channels (Amador et al., 2013, Jiang et al., 2003) (Fig. 1.4) and measurement of single channel currents (Watras et al., 1991) both give evidence sug- gesting that ion channels possess a tetrameric structure. That is, a channel has four subunits, likely identical and each with independent receptor sites. IP3 -sensitive calcium channels have been found experimentally to have an open probability distribution that depends strongly on the present Ca2+ concentration, ex- hibiting a bell-like curve (Fig. 1.5). This implies that calcium has both a stimulating and an inhibitory effect on the channel dynamics: when the calcium concentration at a channel is small, the stimulation is weak and the channel does not open often, whereas when the concentration is large, the inhibitory effect is strong and the chan- nel does not open often. Therefore, it is reasonable to assume that each subunit of a calcium channel has three kinds of receptor sites: sites to which IP3 binds, sites 8 to which Ca2+ binds in such a way so as to encourage activation of the channel, and sites to which Ca2+ binds in such a way so as to inhibit activation of the channel. When a single isolated channel is ac- tive, it opens and closes repeatedly in rapid succession, exhibiting a bursting pattern (Swillens et al., 1998, R¨ udiger et al., 2007, Dupont et al., 2011, Skupin et al., 2010). These bursts are referred to as blips throughout the literature Figure 1.5: The fraction of time a single and typically last around half a second channel is open as a function of local cal- cium concentration. Includes both exper- (R¨ udiger et al., 2007) with individual imental data from (Mak et al., 1998) and simulation results from (R¨ udiger et al., openings within a blip lasting around 2- 2007). Figure recreated from (R¨ udiger 3 ms (Watras et al., 1991, Dupont et al., et al., 2007). 2011, Swillens et al., 1998). Individual channels, however, make up a very small part of the overall picture. While an individual calcium channel has a radius of 6 nm (Skupin et al., 2010, R¨ udiger et al., 2007, Dupont et al., 2011, Swillens et al., 1998), a cell measures about 20 µm (Dupont et al., 2011). Additionally, channels group into clusters around 60 nm across (Swillens et al., 1999), which contain about 10-30 channels (Skupin et al., 2010, Swillens et al., 1999) each. The typical distance between clusters in a cell is around 1-3 µm (Dupont et al., 2011, Swillens et al., 1999). The proximity of channels within a cluster couples their dynamics; one channel opening increases the calcium concentration at nearby channels, which - depending on the exact dynamics - may provoke these other channels to open also. When several channels in one cluster open together in this fashion, it is called a puff. Ex- perimental evidence suggests that puffs typically involve the opening of around five channels (Bootman et al., 1997, Parker and Yao, 1996), and that they actually pre- 9 dominate in clusters (Sun et al., 1998, Thomas et al., 1998). That is, puff events happen far more frequently than blip events - as many as 95% of cluster events are puffs. Swillens et al. (1999) offer compelling arguments based on these results for the internal geometry of clusters. In particular, neighboring channels in a cluster are likely directly adjacent to one another. Thus far, we have laid out the biological details of calcium-induced calcium re- lease. The rest of this chapter will be devoted to discussing the general approach to modeling different aspects of the overall problem. 10 1.4 At the channel: channel stochastics In the previous sections, we de- scribed the internal structure of an ion channel. Specifically, an IP3 -sensitive calcium channel consists of four identi- cal subunits, and each subunit has three different kinds of receptor: IP3 recep- tors, activating Ca2+ receptors, and in- hibiting Ca2+ receptors. In this section, we will describe the overall approach to modeling the chemical dynamics of a sin- gle channel. A single molecule binding or disso- ciating from a receptor site can provoke a conformational change for the entire Figure 1.6: The original De Young-Keizer channel, causing it to open or close. Be- model. The quantities aj are binding cause of this sensitivity, and because of rates and the bj are dissociation rates (Ta- ble 1.1). Figure recreated from (De Young the naturally discrete state space, an ob- and Keizer, 1992). vious way to model channel transitions is as a stochastic process, defining each possible combination of bound receptors as a different state, defining transition rates between them based on the chemical reac- tion rates for the receptor binding and dissociation, and defining an active subunit state and an opening threshold for the channel: when the number of subunits simul- taneously in the active state reaches this threshold, the channel opens and releases calcium. When the number of active subunits falls below this threshold, the channel closes again. 11 a1 400 µM−1 s−1 IP3 binding with no inhibiting Ca2+ bound b1 52 s−1 IP3 dissociation with no inhibiting Ca2+ bound a2 .2 µM−1 s−1 Ca2+ binding to inhibiting site with IP3 bound b2 .2098 s−1 Ca2+ dissociation from inh. site with IP3 bound a3 400 µM−1 s−1 IP3 binding with inhibiting Ca2+ bound b3 377 s−1 IP3 dissociation with inhibiting Ca2+ bound a4 .2 µM−1 s−1 Ca2+ binding to inhibiting site with no IP3 bound b4 .0289 s−1 Ca2+ dissociation from inh. site with no IP3 bound a5 20 µM−1 s−1 Ca2+ binding to activating site b5 1.65 s−1 Ca2+ dissociation from activating site Table 1.1: Channel transition rates from (De Young and Keizer, 1992). Refer to Figure 1.6 for a schematic illustrating the behavior of channel transitions. The simplest and most widely used version of this kind of model, introduced by De Young and Keizer (1992), consists of one of each of the three kinds of receptor on each subunit. A subunit is considered to be in an active conformational state when IP3 is bound and calcium is bound to the activating site but not to the inhibiting site. The channel is considered to be open when at least three of the four subunits are active. In the De Young-Keizer model, the dynamics for binding and dissociation of molecules at the receptor sites are assumed to be noncooperative between subunits. That is, the rate constants for each binding or dissociation event are dependent only on the local calcium and IP3 concentrations, and possibly on the current state of the subunit, but not on the states of the other subunits (Fig. 1.6). Specifically, there are eight states in which a subunit can be, and we index these states with a 1 for each molecule if it is bound and a 0 if it is not, in the order IP3 , activating Ca2+ , inhibiting Ca2+ . For instance, the state S110 is the one in which IP3 is bound and Ca2+ is bound to the activating site but not to the inhibiting site, and S100 is the state in which only IP3 is bound. Statistically, as with any combination reaction, the rate at which binding events occur depends on the number of the binding species present, and so the binding rates 12 for IP3 and Ca2+ depend on the concentrations of these species near the channel re- ceptors. See Table 1.1 for the values of the rate constants De Young and Keizer used in their model. In Chapter 3, we will examine the choices of rate constants used in different models and elaborate on a method, originally developed by Swillens et al. (1999), for determining the rate constants from experimental data. 13 1.5 Between channels: reaction-diffusion The geometry of the cell, and even the internal geometry of a single cluster, is on a much larger scale than the activity of individual calcium ions, and so it is natural to assume a continuum approximation of the calcium ions away from the channels and model the concentration with a standard diffusion equation: ∂c = D∇2 c ∂t where D is the diffusion constant for calcium and c = c(r, t) is the calcium concen- tration at the point r and at time t. The boundary conditions for this equation are chosen depending on the scope of the model; a model encompassing the entire cell will necessarily take into account physical conditions at the cell boundary (often either no-flux conditions or a specific flux corresponding to an extracellular signal), whereas a model focusing on a sin- gle channel might choose a Dirichlet condition at some arbitrary distance from the channel which is much smaller than the cell length. We briefly note that it is not immediately clear that a continuum approximation makes sense: if we consider a box whose side length is the diameter of a channel (Rch = 6 nm), then a concentration of [Ca2+ ] = .05 µM corresponds to an average of 5.2 · 10−5 ions in this box. Additionally, given a channel current of .1 pA and a mean channel open time of 3 ms, every time a channel opens, it releases only about 950 ions. Over the course of a blip event, this number might rise by an order of magnitude. These numbers collectively seem far too small for a continuum approx- imation to work except far from any channel, where the action of individual ions is unimportant. Nevertheless, the majority of work in the literature uses a continuum approximation, even in the region where a channel is situated. 14 The most important addition to the diffusion equation is, of course, the influx of calcium when a channel is open. To model this influx, source terms are added to the equation. There are of course many choices for the form of these source terms, but the most common choice is to model them by point sources: ∂c X = D∇2 c + φj (t)δ(r − rj ) ∂t j where rj is the location of the j th channel, δ is the usual Dirac delta,   φ0 , j th channel is open at time t  φj (t) =  0, j th channel is closed at time t  and φ0 is the single channel flux strength. In Section 2.4, we will examine and discuss different choices for modeling the source terms. Other molecules exist naturally in the cytoplasm that can bind to calcium and prevent it from binding to channel receptor sites. These molecules are commonly referred to as “buffers”, “ligands”, or “chelators” in the literature. The effect of buffers is quantitatively significant for the diffusion of calcium, and so is also typi- cally included in models. Additionally, experiments probing the behavior of calcium channels will often involve adding a dye that binds to calcium ions, allowing the researcher to measure more easily the concentration of calcium in different locations. Thus, in order to be able to compare the result of a model quantitatively to those of an experiment, it is necessary to include the effects of the dye in the model. Naturally-occuring buffers are called endogenous and include mobile buffers which diffuse and stationary buffers which do not. The dye molecules are referred to as dye buffers or exogenous buffers. Inclusion of these buffers in the model results in additional reaction terms in the equation for c as well as introducing reaction-diffusion equations for each species of 15 buffer: ∂c X  X = D∇2 c + ki− (Bi − bi (r, t)) − ki+ bi (r, t)c(r, t) + φj (t)δ(r − rj ) ∂t i=s,m,e j ∂bs = ks− (Bs − bs (r, t)) − ks+ bs (r, t)c(r, t) ∂t ∂bm − = Dm ∇2 bm + km (Bm − bm (r, t)) − km + bm (r, t)c(r, t) ∂t ∂be = De ∇2 be + ke− (Be − be (r, t)) − ke+ be (r, t)c(r, t) ∂t where ki± are the reaction rates of the buffers with calcium; Bi is the total concentra- tion of the ith buffer species; bi is the free (unbound) concentration of the ith buffer species; and the indices s, m, and e refer to the stationary, mobile, and exogenous (dye) buffers, respectively. The total and free concentrations of buffer, B and b, will in general depend on location and time. Generally, the total concentration is much larger than the bound concentration (this also means the free concentration is much larger, and this will help motivate Section 2.3). Additionally, the constant mixing between free and bound buffers means the gradients of total buffer concentrations ought to be much more subtle, even in a region of high activity. Therefore, most models simply take B to be constant. Since the primary effect of buffers is the slowing of transport of calcium away from an open channel, some models aiming for simplicity will eschew explicit repre- sentation of buffers in the model, opting instead for a reduced “effective diffusion” constant for the calcium ions. In Section 2.3, we will examine the strength of the effect of buffers compared to other parts of the model, and explore other simplified methods for the treatment of buffers. Models aiming for thoroughness may include other terms or even whole equations for other effects. For instance Skupin et al. (2010) use a separate variable E for the lumenal calcium concentration (i.e., that in the ER) which also varies in location and time, and let the current of an open channel depend on the value of (E − c) at the 16 channel. Skupin’s model additionally includes terms modeling pumps transporting calcium back into the ER from the cytoplasm. These effects are beyond the scope of this thesis, but are mentioned here for completeness. 17 1.6 Coupling the regimes: the Gillespie algorithm The dynamics of calcium in the cytoplasm and the dynamics of individual molecules and ions binding to receptor sites take place on drastically different spa- tial and temporal scales. Recall (see Section 1.3), the radius of a channel is 6 nm, while the typical distance between clusters is 1-2 µm, and the cell is about 20 µm across. The mean open time for a single channel (which rapidly opens and closes while bursting) is around 3 ms, while the timescales for whole blips or puffs can easily span hundreds of milliseconds (R¨ udiger et al., 2007, Swillens et al., 1999), and larger events (especially spikes) can exhibit behavior on timescales of minutes (Skupin et al., 2010). In contrast, the average time between transitions for a single channel can be as small as 1 ms, and the inclusion of hundreds of channels throughout the cell means transitions throughout the cell are happening on a scale of microsec- onds. Therefore it makes sense to model the large-scale diffusion and the small-scale channel receptor dynamics separately and find a way to couple the models together. While diffusion of calcium in the cell can naturally be modeled in a continuum approximation using partial differential equations, each channel is in one of a finite, discrete set of states, and the transitions between states are much better suited to probabilistic modeling. On the other hand, centering the model on the channel transitions can prove com- putationally expensive. A naive approach might be to pick a timestep sufficiently small that all channel transitions are rare events, and model it as a traditional Markov chain for a Poisson process. The drawback is that computation of diffusion then has to be done at this timestep, which is much smaller than that needed to model the diffusion accurately. This motivated the development of a hybrid scheme by R¨ udiger et al. (2007) to strike a balance between these two disparate regimes. The goal is to maximize 18 the amount of simulation time between stochastic computations so that the timestep used to compute diffusion, which is necessarily bounded by the time between stochas- tic events, may be relaxed. The technique relies on the adapation of an algorithm developed by Gillespie (1977) which consists of drawing two random numbers to determine the properties of the next transition: let αj denote the transition rates P away from the current state of the system, and let α ≡ j αj . Draw two random numbers r1 , r2 ∈ U (0, 1). Then the time to the next transition, τ , is given by   1 1 τ = log α r1 and the index j of the next transition is determined by j−1 j 1X 1X αk < r2 ≤ αk α k=1 α k=1 See Appendix A for a derivation of these formulas. One will notice that this algorithm creates nonuniform intervals between channel transitions, but this typical transition time is much larger than the timestep of the naive Markov-Poisson model, because each timestep is guaranteed to see precisely one transition instead of having a transition being a rare event. This means that the upper bound the Markov model imposes on the diffusion timestep is much higher, and so a larger diffusion timestep may be chosen with efficiency in mind. The adaptation by R¨ udiger focuses on the fact that, because the channel transi- tions depend on the local calcium concentration at the channel’s receptor sites, and because the calcium concentration is necessarily often in flux, the transition rates themselves are not constant in time. In other words, Gillespie’s algorithm assumes that all of the transition rates are constants, while R¨ udiger’s implementation ex- tends the algorithm to time-dependent rates as follows: let αj = αj (t, c(t)) and let α = α(t, c(t)) be defined analogously as before. Then the time to next transition is 19 given by (see Appendix A for details) ˆ τ   0 0 0 1 α(t , c(t ))dt = log 0 r1 Note that, in the event that α is constant, this equation reduces precisely to   ´t τ = α1 log r11 . If we define g(t) = 0 α(t0 , c(t0 ))dt0 , this suggests the following initial value problem: dg = α(t, c(t)); g(0) = 0 dt   which can be evolved numerically until g reaches the value log r11 . Since αj > 0,   a solution g(τ ) = log r11 is guaranteed to be unique. Existence of a solution is not mathematically guaranteed for every set of αj , but such situations are decidedly un- physical. In particular, in the context of modeling calcium dynamics, such a situation would require all of the IP3 and Ca2+ to leave the cell cytoplasm. Thus, Gillespie’s algorithm may still be applied, in R¨ udiger’s form, to the problem of modeling cal- cium ion channel dynamics, allowing the optimization of the stochastic timestep. 20 1.7 Summary In this chapter, we have outlined the biological phenomenon of calcium-induced calcium release as well as the general approaches to modeling diffusion of calcium away from channels and modeling channel receptor reactions. We have also discussed a method for coupling together these two regimes in a precise way. In the next two chapters, we will delve more deeply into the details of employing these modeling approaches and examine different choices that may be made, along with their advantages and drawbacks. We will also examine some approximations that can be made to simplify the model and reduce computation time, and discuss the relative impact of these approximations on the accuracy of the overall model. The purpose of our work was to streamline existing approaches to modeling CICR, that is, to improve the overall accuracy and efficiency by examining different modeling approaches. The specific results we found are as follows: • The assumption that the calcium concentration immediately equilibrates to a steady-state spatial profile upon the opening or closing of a single channel is a useful and reasonable approximation. • The assumption that the concentrations of any buffer species do not depend on location or time is a useful and reasonable approximation. • It is possible to model CICR treating a cluster as a single, fundamental stochas- tic object and obtain qualitatively reasonable results. • The assumption that the single-channel flux is concentrated at a point is not a reasonable approximation and is of questionable usefulness. Instead, taking the flux to be uniform over the channel volume offers an approach with less 21 arbitrariness that is neither computationally expensive nor difficult to imple- ment. • It is possible to construct channel receptor binding and dissociation rates for a cluster of channels using data from experimental results. The use of such a construction yields a working numerical model for a cluster of calcium channels. The first three results are significant for their implications towards model effi- ciency; these simplifications dramatically reduce computation time and at least the first two do so without significantly impacting the accuracy of the model. The fourth result is important due to its ubiquity in the literature. We examine the drawbacks to this approach and offer an easy alternative. Finally, the last result is, to our knowledge, a novel result that may allow for more careful analysis of the inner work- ings of calcium channel binding reactions as well as ultimately provide the means for more exact quantitative models of CICR. Chapter TWO Approximations 23 2.1 Overview With any computation of sufficient complexity, a choice must be made in how to strike a balance between thoroughness and computational tractability. Many mod- els in the literature focus on examining a single cluster of calcium channels or just one isolated channel. These problems involve sufficiently few quantities that a fully accurate model is not very computationally expensive. However, when one considers an entire cell with as many as 50 clusters, the cost goes up dramatically: diffusion occurs on a much wider domain but must still be treated with sufficiently fine granu- larity to account for the internal geometry of a cluster, the spherical symmetry used to model diffusion is lost, the number of possible transitions increases tremendously, and the expected time between channel transitions becomes one or two orders of magnitude smaller. Computational techniques exist to mitigate these problems to a degree. Never- theless, models examining an entire cell typically require the use of a supercomputer to become manageable. One of the primary aims of this thesis is to examine approx- imations and simplifications that can dramatically reduce the computational time for a full-cell model while keeping sufficient accuracy to provide at least qualitatively useful results. In this chapter, we will discuss timescales of different parts of the model, the relative impact of buffer dynamics, several simplifications involving the limiting so- lution of an open channel, and the relationship between single channel statistics and cluster statistics. In particular, we will examine the following approximations: • the assumption that the effect of channel influx occurs instantaneously • the assumption that the buffer concentrations are independent of location and time 24 • the assumption that the flux due to an open channel is concentrated at a point • a simplification wherein puffs, instead of blips, are the fundamental stochastic event and individual channel behavior is not explicitly computed We will demonstrate that the first two assumptions are good approximations, that the third is not, and that the fourth is plausible. 25 2.2 Instantaneous equilibration 2.2.1 Motivation When recreating existing models, it quickly became apparent that the computa- tions required to calculate the diffusion of calcium ions dwarfed those for the channel transitions in terms of number of computations per timestep. This proved to be true regardless of the approach used in computing the diffusion. Skupin et al. (2010) ap- proached the problem by using Green’s functions in order to be able to compute the calcium concentrations only locally at each channel location instead of everywhere on a large mesh. However, because of the geometry of the problem, the Green’s func- tions involved infinite series of eigenfunctions, the approximation of which still left many terms to be computed and proved to be computationally expensive. Skupin et al. made explicit mention of the use of a supercomputer to analyze their model. It seemed apparent that the most dramatic way to speed up the computation time was to find a way to eschew the diffusion computation altogether. When modeling a single calcium channel, we noticed that the calcium concentration increased rapidly near a channel upon opening and quickly leveled off. Additionally, this profile had limited spatial extent; beyond a certain distance - only a fraction of the cell radius - the difference between c(r) and the base cytosolic concentration c0 was a very small fraction of c0 . In models with relatively high buffer concentrations, we even found that every- where in the cell, the spatial calcium profile approached a steady state, and that after a fraction of the cell radius, the spatial profile dropped to within machine error of c0 . A single channel opening or closing therefore had a very rapid effect on an in- tracluster scale and a very weak effect on an intercluster scale. Thus, we examined 26 the possibility of assuming that a channel instantaneously achieves this steady-state profile upon opening, and that the profile instantaneously decays back to c0 upon closure of the channel. 27 Figure 2.1: Time evolution of the calcium concentration at an open channel with no approximations. Biological parameters used are from (R¨ udiger et al., 2007). Modeled using a finite-difference scheme with the full reaction-diffusion equations for c and bi with three species of buffer. Channel flux distribution uniform over the channel volume. 2.2.2 Fast timescale of channel influx The three possible mechanisms for closing a channel are: 1. dissociation of IP3 from an active subunit, 2. dissociation of Ca2+ from the activating site of an active subunit, and 3. binding of Ca2+ to the inhibiting site of an active subunit. Since three active subunits are required for a channel to be open, the average closing rate will then be roughly three times the sum of the rates for these three events. With the nomenclature of Table 1.1, the closing rate will be 3 · (b1 + b5 + a2 c) and the average time to close will be the reciprocal of this quantity. Looking at the parameters from Table 2.1, and taking an at-channel peak concentration of 50 µM, the fastest average 28 1 closing time possible is roughly tclose = 500 s−1 = 2 ms, which also matches well the observed mean open time for a single channel (Section 1.3). Thus, if the concentration at the channel levels off or even slows substantially in much less than a millisecond, a channel that has just opened is extremely unlikely to undergo a state transition before the local concentration reaches this plateau. Thus for the purposes of computing channel transitions, it would be a good approximation to assume the concentration reaches this level instantaneously, as doing so would make no difference to the channel’s dynamics. We found such a rapid leveling off to be the case for a variety of sets of biological parameters and stochastic transition rates. In Figure 2.1, we have shown this leveling off for a model consisting of a single channel in the center of a cell- sized sphere with Dirichlet boundaries, giving a full dynamical treatment of the buffers. It is evident that in this case, the concentration doesn’t quite reach an Figure 2.2: Time evolution of the spatial equilibrium value within the first mil- calcium concentration profile. Recreated from (Smith et al., 1996). Of particular lisecond. Nevertheless, it increases ex- note is the 10 mM exogenous buffer con- centration. tremely rapidly for the first 50 microsec- onds or so, and then for the rest of that first millisecond, rises extremely slowly. In situations where the model includes a large buffer concentration or ignores the dynamics of the buffers and takes them to be constant (an approximation that will be explored in depth in Section 2.3), the sit- uation is even more dramatic. In both cases, the spatial concentration profile quickly reaches a steady-state distribution. In the former case (Fig. 2.2), the entire profile equilibrates within 1 µs. In the latter (Fig. 2.3), the entire profile equilibrates within .6 ms, and the at-channel concentration equilibrates within .1 ms. In any case, the 29 Figure 2.3: Time for the concentration at a point to reach 99.9% of the equilibrium concentration at that point. All model parameters are the same as in Figure 2.1 ex- cept the buffer species concentrations are treated as constant - both uniform through- out the cell and independent of time. The reason the time falls to 0 at around 1.3 µm is because c0 is already within 99.9% of the equilibrium value. Beyond this point, the effect of the channel is negligible. effect of an open channel on itself can justifiably be assumed to be instantaneous with respect to the receptor binding reactions. 30 (a) (b) (c) (d) (e) a1 20 400 80 20 ki+ 35.8 b1 20 52 .64 52 ki− 5 a2 .001 .2 .04 .001 kd+ .088 b2 .03 .2098 .48 3.077·10−4 kd− 1.67 a3 2.6 400 80 20 ki+ 35.8 b3 20 377 40 52 ki− 5 a4 .025 .2 .4 .025 kd+ .088 b4 .1 .0289 .0768 .0025 kd− 1.67 a5 10 20 15 10 ka+ 166 b5 1.225 1.647 12 122.5 ka− 65.9 Table 2.1: Comparison of channel rate constants from various models. Units are µM−1 s−1 for binding rates (ai ) and s−1 for dissociation rates (bi ). Repeated and elaborated upon in Chapter 3 (Table 3.2). 31 Figure 2.4: Spatial distribution of calcium due to a cluster with five open channels 2.5 ms after opening. 2.2.3 Limited spatial extent of channel influx In the previous section, we illustrated that a channel opening has an extremely rapid effect in its vicinity. We also touched on the fact that, although the effect of an open channel is felt more slowly far from the channel, that effect is much smaller. In this section, we will briefly illustrate this fact. The typical distance between clusters is 1-2 µm. We modeled the effect of a source five times as strong as an individual channel (Fig. 2.4) and found that even after a couple milliseconds had elapsed, the effect of the active cluster had only just reached 1 µm. Even at half this distance, the elevated concentration was only about three times the base level - a fraction of the effect of a channel on its neighbors. In fact, if we include the approximation of treating buffers as constant, the spa- 32 tial profile achieves within a couple milliseconds a steady-state profile that does not change to within machine error. This profile drops to within 0.1% of c0 at 1.3 µm from the cluster (Fig. 2.3), and to within machine error of c0 at about 4 µm from the cluster - only 20% of the cell diameter. In the previous section, we demonstrated that the effect of an open channel on itself and its neighbors propagates on a timescale much shorter than the typical stochastic transition time, and so near an open channel, the calcium concentration profile may be assumed to take effect instantaneously. Here, we have shown that the effect of an open channel outside of its home cluster is so small that making this assumption will not have a profound effect on the dynamics of other clusters. We have therefore justified that taking this approximation for a full-cell model is reasonable. 33 2.2.4 Eliminating time evolution from the diffusion compu- tation We have just demonstrated that assuming instantaneous equilibration of the concentration profile upon a channel opening is a reasonable approximation to make. The obvious implication is that continuous time evolution is no longer needed for modeling the diffusion of calcium through the cell. In this section, we illustrate how this leads to a very simple computation of the calcium concentration at any point in the cell. Recall the reaction-diffusion equations for c and b: ∂c X  X = D∇2 c + ki− (Bi − bi (r, t)) − ki+ bi (r, t)c(r, t) + φj (r, t) ∂t i=s,m,e j ∂bs = ks− (Bs − bs (r, t)) − ks+ bs (r, t)c(r, t) ∂t ∂bm − = Dm ∇2 bm + km (Bm − bm (r, t)) − km + bm (r, t)c(r, t) ∂t ∂be = De ∇2 be + ke− (Be − be (r, t)) − ke+ be (r, t)c(r, t) ∂t The first equation is affine in c and we can exploit this fact in combination with instantaneous equilibration to make computation of the calcium concentration at any point in the cell quite easy. Suppose we have two steady-state solutions c1,2 (r) corresponding to open channels at two distinct locations r1,2 . For simplicity, assume only one buffer species. Then cj satisfies: D∇2 cj − k + b(r)cj (r) + k − (B − b(r)) + φ(r − rj ) = 0 34 This equation has a constant solution c0 = K B−b b0 0 when all channels are closed (see Appendix B for more details). Let cj (r) ≡ c0 + fj (r). Then the equation reduces to D∇2 fj − k + b(r)fj (r) + φ(r − rj ) = 0 So f ≡ f1 + f2 satisfies D∇2 f − k + b(r)f (r) + φ(r − r1 ) + φ(r − r2 ) = 0 and so c ≡ c0 + f1 + f2 = c0 + (c1 (r) − c0 ) + (c2 (r) − c0 ) satisfies D∇2 c − k + b(r)c(r) + k − (B − b(r)) + φ(r − r1 ) + φ(r − r2 ) = 0 which is the full time-independent equation for two open channels. There is a subtlety here, which is that c1 and c2 actually satisfy equations with different functions b(r) (recall, however, that we assume from the start that B is con- stant). Neher (1986) demonstrates that the maximal change in buffer concentration due to an open channel is quite small compared to the total buffer concentration, and so takes the buffer concentration to be uniform (an approximation that will be looked at more closely in the following section). This assumption implies the above result. ∂cj If c1 and c2 satisfy either cj = c0 or ∂n at the cell boundary (which are the two most realistic scenarios), then clearly so will c, and so c is also a solution of the full boundary-value problem. More generally, this means that if ck (r) is the steady-state solution for the k th channel (which may, in general, depend on the geometry of the cell), then the full 35 steady-state solution for an arbitrary number of channels N is given by N X c(r) = c0 + (ck (r) − c0 ) · ok k=1 where   1, k th channel is open  ok =  0, k th channel is closed  which is trivial to implement numerically. If we assume instantaneous equilibration as detailed earlier in this chapter, then we can allow the ok to depend on time, and the full time-dependent solution for the problem becomes N X c(r, t) = c0 + (ck (r) − c0 ) · ok (t) (2.2.1) k=1 In other words, assuming instantaneous equilibration allows us to reduce solving the PDE to what amounts to a single inner product of two vectors, one of which does not depend on time and can be precomputed solely from considerations of geometry. 36 2.2.5 Handling the geometry With the use of equation 2.2.1, the time evolution of diffusion is completely removed from the computational process, and the only aspect of solving the reaction- diffusion equation that remains is examining the interaction between the sources and the cell boundary. We examine a few easy choices. If we assume Dirichlet conditions at the cell boundary, then this formula is very easy to implement generally over the whole domain of the cell. A spherical domain with a constant value on the boundary has a simple solution using the method of images: for each channel in the cell, we can create a corresponding image source outside the cell with an appropriately scaled strength. The drawback to this method is that a Dirichlet condition on the cell boundary is unphysical and likely to be a poor approximation to a real biological situation. A Neumann boundary condition is much more physically appropriate, and can still be treated using the method of images, but spherical symmetry doesn’t allow for such a simple, elegant solution to this problem. Instead, the solution yields a Green’s function consisting of an infinite series of terms, much like that used by Skupin et al. (2010), and so any computational efficiency gained from employing the instantaneous equilibration assumption is lost again in approximating this series. As demonstrated in Section 2.2.3, the single-channel steady-state profile has a finite extent depending on the biological parameters chosen for the model; using the set of biological parameters from (R¨ udiger et al., 2007) and treating the buffer concentrations as constant, we had c(r) − c0 < mach for r > .35 Rcell . This lends itself to another possible approximation for the full cell geometry, namely restricting the possible locations of the channels to a region in the center of the cell this distance from the cell boundary. This situation is further justified by the possibility of IP3 - sensitive channels being localized in a region roughly 6 µm wide (Callamaras and 37 Parker, 1999) while the cell is roughly 20 µm across. In this case, the effects of any open channel therefore never make it as far as the cell boundary, and we can use the simplified formula for c(r, t) everywhere in the cell without worrying about the boundary condition; both a constant boundary condition and a no-flux boundary condition are automatically satisfied to within machine error. In this case, if we take instantaneous equilibration of channel dynamics as an assumption, then the full solution to the reaction-diffusion equation for c becomes N X c(r, t) = c0 + c(r − rk ) − c0 ) · ok (t) (˜ k=1 where c˜ is the steady-state solution for a single channel at the origin and rk is the location of the k th channel. That is, only one steady-state distribution needs to be computed in this case, along with the set of pairwise distances between channels. This is clearly a remarkably simple and computationally inexpensive method to im- plement in a numerical model. 38 2.2.6 The Gillespie-R¨ udiger algorithm In the last section, we described how the assumption that an open channel in- stantaneously reaches a steady-state spatial distribution, along with the assumption of either of two realistic boundary conditions, leads to a remarkable simplification of the computation of diffusion of calcium throughout the cell. While this is likely the most dramatic impact of this approximation - because the computation of diffusion is the most computationally expensive part of the model - it is not the only one. Recall the Gillespie-R¨ udiger algorithm for computing the stochastic transitions: P given stochastic transition rates αk (t) with α(t) ≡ k αk (t), we may compute the time of the next transition, τ , by evolving dg = α(t, c(t)); g(0) = 0 dt 1  until g(τ ) = log r where r is a random number selected uniformly from (0, 1). (See Section 1.6 for an overview of the algorithm or Appendix A for details.) Following this algorithm therefore requires numerically evolving an ODE at every stochastic timestep. However, if we use the approximation of instantaneous equili- bration, the transition rates αk no longer depend on time. Between transitions, they are all constant since c is independent of time. Therefore, the computation of τ reduces to the original version given by Gillespie:   1 1 τ = log α r Thus, the numerical evolution of an ODE is replaced by a very small number of operations. 39 2.3 Constant buffers 2.3.1 Motivation Recall the full reaction diffusion system: ∂c X  X = D∇2 c + ki− (Bi − bi (r, t)) − ki+ bi (r, t)c(r, t) + φj (r, t) ∂t i=s,m,e j ∂bs = ks− (Bs − bs (r, t)) − ks+ bs (r, t)c(r, t) ∂t ∂bm − = Dm ∇2 bm + km (Bm − bm (r, t)) − km + bm (r, t)c(r, t) ∂t ∂be = De ∇2 be + ke− (Be − be (r, t)) − ke+ be (r, t)c(r, t) ∂t Here, there are four (or more, if more than one species of each kind of buffer is modeled) separate quantities to compute, and the computational cost of diffusion - which is already the most expensive part of modeling CICR - scales with the number of these quantities. Being able to reduce this system to one equation would therefore have a positive impact on the computational cost of the model, and because the calcium concentration is the principle quantity of interest, the simplest way of doing this is simply to ignore the other three equations and model the unbound buffer concentrations as uniform throughout the cell and constant in time (recall that the total buffer concentrations are already assumed to be thus). Additionally, as alluded to in the previous section, this assumption of constant buffers couples nicely with the assumption of instantaneous equilibrium, making the implementation and justification of the latter significantly easier. In particular, with the assumption of constant buffers, the time evolution of the spatial calcium pro- file due to influx of an open channel not only slows rapidly, but in fact reaches a steady-state distribution on a timescale that is on the same order of magnitude as 40 the average stochastic transition time for a single channel. Using the assumption of constant buffers allows us to model a single channel with any set of parameters and any source distribution desired, quickly compute the steady-state distribution, and use that distribution throughout the rest of the model. The primary impact of buffers is twofold: altering the diffusion dynamics of the calcium concentration, and affecting the equilibration speed at an open channel. In the next section, we will demonstrate that the assumption of constant buffers does not significantly influence either of these effects, and that it is therefore a reasonable assumption to make. 41 Figure 2.5: Comparison of the spatial profiles due to the influx of a single channel with (a) constant and (b) dynamically evolved buffers. Elapsed time is 2.5 ms. Biological parameters from R¨ udiger et al. (2007), with uniform flux over channel volume. Computed with 2nd-order radially symmetric finite difference scheme and 4th-order Runge-Kutta method. Dirichlet conditions at cell boundary. See Figure 2.6 for a plot of the absolute error. 2.3.2 Justification The primary effect of buffers is on the diffusion dynamics of calcium. In Sec- tions 2.2.2 and 2.2.3, we illustrated that there are effectively two regimes for diffu- sion: far from a channel, where the effect of influx due to an open channel is felt slowly but weakly, and near the channel, where this effect is felt strongly but rapidly. In this section, we will demonstrate that the assumption of constant buffers preserves the weakness of the former and rapidity of the latter. Figures 2.5 and 2.6 compare the results of modeling a single channel with con- stant or dynamic buffers. The elapsed time in both cases is 2.5 ms, which (see 42 Figure 2.6: The absolute error corresponding to Figure 2.5. The error is small compared to c0 = .05 µM far from the channel. Section 1.3) is roughly equal to the expected open time for the channel. Therefore, even though, should the channel remain open, the dynamic-buffer profile would con- tinue to rise (Fig. 2.1), the channel is unlikely to remain open much longer than this amount of time. These results show that far from an open channel, the effect of the constant buffer approximation is negligible and preserves the weakness of the open channel’s effect at these distances. Next, we will examine the effect of this assumption at the channel itself. Sup- pose that the buffer concentrations are indeed constant, and for simplicity, assume only one buffer species. Additionally, consider a situation with only one channel. The reaction diffusion equation for calcium is then ∂c = D∇2 c − k + bc(r) + k − (B − b) + φ(r − rj ) ∂t 43 When the channel is closed, there is a constant solution c0 = K B−b b , and when the channel is open, assume there is a steady-state solution c(r) = c0 + c˜(r) with c˜ satisfying D∇2 c˜ − k + b˜ c(r) + φ(r − rj ) = 0 The steady-state solution in the case that φ is a point source is derived in Appendix B. Then the full time-dependent solution when the channel is closed is given by cclosed = c0 + F (r, t) and the solution when the channel is open is given by copen = c0 + c˜(r) + G(r, t) where F and G both satisfy ∂F = D∇2 F − k + bF ∂t This equation is separable: let F (r, t) = R(r)T (t). Then we have: 1 dT 1 D∇2 R − k + bR ≡ −ν  = T dt R The equation for T yields T (t) = e−νt and so ν is going to be a characteristic decay rate for a particular mode. The modes are determined by D∇2 R + (ν − k + b)R = 0 along with boundary or regularity conditions. This equation will only have solutions that are regular over the whole domain when ν −k + b ≥ 0 and so the smallest possible 44 (a) Stationary buffer (b) Mobile buffer (c) Exogenous buffer Figure 2.7: The time evolution of the buffer species concentrations at an open chan- nel. Total elapsed time 2.5 ms. value of ν is given by ν0 = k + b. Since ν represents a decay rate, the smallest value will represent the slowest decay mode, and thus represent the overall equilibration rate 1 of the solution. In particular, it should be the case that ν0 gives the characteristic equilibration time. Generalizing to several buffer species, it is easy to see that X ν= kj+ bj j 45 and so we may verify this result using, for instance, the biological parameters from (R¨ udiger et al., 2007). In this case, we have bs = 78 µM, bm = 225 µM, be = 39 µM, ks+ = 50 µM−1 s−1 , km + = 5 µM−1 s−1 , and ke+ = 150 µM−1 s−1 . Therefore, ν0 = 11000 s−1 and 1 ν0 = 92 µs, which matches closely the results seen in Figure 2.1. This computation has been performed assuming that the buffers are constant. Removing this assumption, an analytic examination of the decay modes is not so easy, but it should be the case that in a small time interval, the decay rate in that interval is approximately what it would be retaining this formula for ν0 . Allowing the concentrations to evolve numerically for 2.5 ms, using the same model as in Figure 2.5 (with dynamic buffers), we find (Fig. 2.7) the buffer concen- trations to be bs = .6 µM, bm = 224 µM, and be = 17 µM after this amount of time. 1 Using these values, we find ν0 to be 270 µs. This is three times as large as the value obtained using constant buffers, but is still small compared to the mean open time of 2 ms. The assumption of constant buffers, therefore, not only gives a result consistent with the numerical results from a simulation with no approximations (Fig. 2.1), but gives a result that is robust upon removal of the assumption. We have therefore shown that assuming the buffer concentrations are constant is a good approximation both near an open channel and far away from it. 46 2.4 Channel flux distribution 2.4.1 Motivation Because of the self-affecting nature of calcium ion channels, it is important to be able to accurately model the effect of an open channel on itself. Furthermore, the influx from an open channel creates a strong concentration gradient (Fig. 2.5), and so the effect of an open channel on itself is much greater than the effect even on directly adjacent channels. Therefore, it is critical to be able to model the concentration at an open channel accurately. Much of the literature (e.g., Skupin et al. (2010), Swillens et al. (1998)) chooses to model the influx of calcium using a point source. That is, φ(r) = φ0 · δ 3 (r) where φ0 is a constant characterizing the flux strength and δ 3 (r) is the usual three- dimensional Dirac delta. In this section we will examine this choice of channel flux distribution, discuss its implications and drawbacks, and discuss an alternative: the uniform flux  φ0 , r < Rch   Vch φ(r) =  0,  r > Rch with Vch = 43 πRch 3 the volume of the channel and φ0 the same value as for the point source. 47 2.4.2 Localized flux The most common and, perhaps, simplest distribution for modeling the calcium flux due to an open channel is a point source: φ(r) = φ0 · δ 3 (r) The primary advantage to this distribution is that the reaction-diffusion equation ∂c = D∇2 c − k + bc + k − (B − b) + φ0 · δ 3 (r) ∂t for the calcium concentration due to a single open channel with this flux distribution, located in the center of a large domain, has the analytic steady-state solution q + − kD b r φ0 e c(r) = c0 + 2πD r k− B−b where c0 = k+ b = K B−b b . A derivation of this solution may be found in Ap- pendix B. This solution, while very easy to implement outside a channel, has the drawback of becoming singular at r = 0 and therefore has a wide (indeed, infinite) range of values within the domain of the channel. However, we know that the dynamics of the channel depend strongly on the local calcium concentration (Fig. 1.5). Conse- quently, it is important to find a method to select the value of the concentration at the channel. Skupin et al. (2010) do not use the analytic solution, but nevertheless assume lo- calized sources for the purpose of exploiting their ease of use with a Green’s function approach to modeling the diffusion of calcium and buffers through the cytoplasm. As a result, the singularity at the channel is something they were forced to address 48 Figure 2.8: Steady state-concentration due to influx of a single open channel with the flux (a) spread uniformly over the channel volume, (b) localized at the channel center. For the localized source, the concentration profile is truncated at the channel boundary. in the implementation of their model, and they chose to assign the value of the con- centration at an open channel equal to the value at a distance some fraction of the channel radius away from the center of the channel. This approach is rather difficult to use, as one must justify the distance cho- sen: a typical concentration at the boundary of an open channel might be 30 µM (Fig. 2.8). Comparing this value to Figure 1.5, we see that even with this range of possible values, there is wide variation in the possible dynamics, and therefore a value must be selected without arbitrariness, and this is impossible to accomplish solely considering localized sources. Thus, although a point source seems at first glance to be mathematically simple to use, it proves to be a difficult choice when accuracy is taken into consideration. 49 Figure 2.9: Steady state profile due to an open channel with uniform flux using a spatial granularity of (a) Rch , (b) R2ch , (c) R5ch , (d) R10ch . 2.4.3 Uniform flux Due to the difficulties in implementing a point source for the channel flux, we therefore examined spreading out flux uniformly over the domain of the channel. For simplicity (in order to maintain spherical symmetry), we chose to let the flux be nonzero over the region r < Rc where Rc is the channel radius. That is, we replaced the flux density term φ0 δ 3 (r) with one of the form   Aφ0 , r < Rch  φ(r) =  0, r > Rch  50 To find the value of A, we recognize that the total flux over the volume should be the same in both cases, and so we easily find that A = (Vch )−1 where Vch = 34 πRch 3 . This is very straightforward to implement numerically, and, with a sufficiently fine spatial granularity, reproduces the analytic profile outside of the channel with a great deal of accuracy (Fig. 2.8). Inside the channel, the profile remains finite and smooth with a peak at r = 0 that we can then select as the concentration at the channel. We briefly note that our choice of flux distribution may be viewed as one of a fam- ily of functions approximating the Dirac delta distribution - that is, it is one member of a Dirac sequence - with compact support on the domain of the channel. One could use the framework of Dirac sequences to systematically study different flux distribu- tions by, for example, constructing a generating function with compact support on the channel domain based on biological considerations (for instance, nonincreasing from channel center to channel boundary), and varying the sequence parameter to match experimental results. See, for instance, (Dang and Ehrhardt, 2012) for a short treatise on Dirac sequences and their generators. It may be tempting to simplify the uniform flux approach back in the direction of a point source by choosing the spatial granularity coarsely enough such that the channel volume comprises the entire central volume element, so that all of the flux occurs in this one volume element. Figure 2.9 shows that this approach (i.e., setting dr = Rch ) introduces up to 30% error. In contrast, the spatial profile has a high Rch degree of accuracy even with a mesh not much finer than this (dr = 5 ). This allows a spatial profile to be computed quite quickly and, when used in combination with the instantaneous equilibration approximation examined in Section 2.2, yields an effective modeling approach for the overall problem. 51 2.5 Puffs as the fundamental event 2.5.1 Motivation Channels in the same cluster are packed closely to one another and so feel the effects of each other quite strongly. Indeed, because of the bursting nature of blips and the short time in which an open channel’s concentration profile decays, it is dif- ficult to see how a puff event could occur were the channels in a cluster not strongly coupled. In fact, experimental evidence (Watras et al., 1991) shows that puffs are the predominant event type in clusters, providing further evidence for this coupling. Although channels in a cluster are strongly coupled, the coupling between clus- ters, is comparatively weak: the typical distance between clusters is around 1-3 µm (Dupont et al., 2011, Swillens et al., 1999), which suggests that the transition from puff level to spike level exhibits a much smaller degree of cooperativity than the transition from blip level to puff level. Because of the strongly concerted action of channels with a cluster, and the stark contrast of this behavior to that at the cellular level, it seems reasonable to treat an entire cluster as one object. Furthermore, a typical cell has a few tens of clusters, a typical cluster has around 10-30 channels, and each channel has four subunits that at any timestep has three possible transitions it can take (or more, depending on the exact model, e.g., (Swil- lens et al., 1999), which includes additional binding sites and cooperativity). This results in an overwhelmingly large state space for the stochastic process. Addition- ally, because a channel’s state depends on the combined states of its subunits, the majority of state transitions don’t cause a channel to open or close, and so do not actually effect a change in the dynamics of the cell. Between the large state space and the strong coupling between nearby channels, it makes sense to cultivate statistics for puffs from a single cluster, and use those 52 statistics to drive cluster dynamics when modeling a full cell. That is, we attempt to construct a full cellular model where each cluster is an isolated, holistic stochastic object, and individual channels are not explicitly included. To do this, we assign to each cluster states corresponding to the number of open channels, ranging from 0 to the number of channels per cluster (typically 15-25). In numerical computations, the state number typically does not exceed 7 because the probability of additional channels opening simultaneously becomes exceedingly rare; the average puff consists of only five channels acting in concert. This reduces the number of states per cluster from around 25 · 4 · 8 = 800 (depending on the model) down to only 8, resulting in a model that is significantly faster. 53 2.5.2 Computation of the cluster transition rates To compute the cluster transition statistics, we model a cluster of individual channels for a sufficiently long simulation time and record the number of open chan- nels at each timestep. From this data, we construct transition rates for the cluster (i.e., rate constants to go from j open channels to k open channels) as follows: Suppose a stochastic process is in a given state, j, out of n total states, and the rate constant to transition from the j th to the k th state is αjk , then the expected time to leave the current state is n !−1 X τj ≡ αjk k=1 and the probability of transitioning to the k th state instead of another is αjk Pjk = Pn = αjk τj k=1 αjk From the data we collect, it is straightforward to extract the time the cluster spends with a certain number of channels open and whether it then opens a new channel or closes one that is firing. It is then simple to compute the average time the cluster spends with j channels firing before opening or closing one - which is precisely τj - and the relative frequency with which it opens a new channel instead of closing a firing one: if nj,k is the number of times the cluster goes from j open channels to k open channels (where k is necessarily either j − 1 or j + 1), then the frequency for opening a new channel is nj,j+1 fj,j+1 ≡ nj,j−1 + nj,j+1 54 and similarly, the frequency for closing an active channel is nj,j−1 fj,j−1 ≡ nj,j−1 + nj,j+1 It is reasonable to assume that this process is ergodic; every state is clearly recurrent, and although under the Gillespie algorithm, the states are periodic (with period 2) by construction, were we to choose a uniform timestep for the Markov chain, there would be a nonzero probability of remaining in the same state (either having not changed, or having gone to another state and back in the same timestep, depending on the size of the timestep), and so the states would be aperiodic. Under this assumption, we have fj,j±1 ≈ Pj,j±1 with equality holding in the limit as simulation runtime goes to infinity. We can therefore derive cluster transition rates by fj,j±1 αj,j±1 = τj provided we run the simulation for a sufficient amount of time. Now that we have computed the cluster transition rates, we can apply the same overall approach we used for modeling the stochastics of channels to the reduced problem where each source is assumed to be an entire cluster. 55 (a) (b) Figure 2.10: A comparison of the resulting dynamics modeling a cluster as (a) con- sisting of 25 channels and using the traditional De Young-Keizer model, (b) a single stochastic object with transition rates as computed in Section 2.5.2. Note: the steep diagonal lines in part (b) result from automatic linear interpolation between two distant data points; the concentration is actually c0 on these intervals. 2.5.3 Results The most straightforward way to verify this approach is to run a simulation with this approach and compare the results to those of a simulation incorporating stochastics of individual channels. Such a comparison may be found in Figure 2.10. Both of these models assume a single cluster in the center of a spherical cell-sized domain with Dirichlet boundary conditions. The biological parameters used are from (Swillens et al., 1999). The states were evolved using the Gillespie algorithm and the calcium concentration was evolved assuming instantaneous equilibrium, and so the total number of open channels was tracked as a function of time (since the spatial extent of the equilibrium profile is much smaller than the radius of the cell, and therefore the average calcium concentration in the cell is precisely proportional to the number of open channels). It is apparent that the results of these two simulations are not dissimilar; the most significant difference between the two is the frequency of puff events, which 56 Figure 2.11: Steady-state spatial profile due to a 5-channel-strength source (typical puff size) with the flux spread uniformly over (a) the volume of the cluster and (b) the volume of a single channel. The radius of a cluster is 30 nm. are roughly twice as frequent in the simplified cluster model (the reason for this discrepancy is still unclear). However, these two schemes show excellent agreement in terms of amplitude and duration of puffs, as well as ratio of puffs to blips. The primary computational expense of the general model, however, is not the stochastic part of the simulation, but the diffusion, if we do not use the assumption of instantaneous equilibration. By averaging the stochastic behavior of a cluster over its entire domain, we treat the cluster as a single source, the strength of which depends on the stochastic state of the cluster. The radius of a cluster is an order of magnitude larger than the radius of a single channel, and we can exploit that by averaging the cluster source over its volume - instead of over a volume the size of a channel - and then choose a much coarser mesh for modeling the diffusion. Figure 2.11 shows that, outside the cluster radius, the two spatial profiles agree closely. 57 Figure 2.12: A spike. Finally, we employed this approach to modeling a cluster in a simulation of a cell with 50 clusters in it. Notably, this simulation - modeling an entire cell - was run in less than an hour on a laptop. An excerpt of the resulting cellular concentration profile as a function of time may be seen in Figure 2.12. Because the average puff event peaks at approximately five channels (which is consistent with Figure 2.12), the one event with 15 open channels is likely due to the simultaneous activation of three separate clusters. In other words, this event is likely a spike. Although this simulation saw only one spike in over two hours of simulation time (compared the expected timescale of minutes), it provides a proof of concept for this approach and demonstrates that CICR spikes may be simulated with a computationally inexpensive approach. Future work may adapt this sufficiently to produce qualitatively or even quantitatively correct results. Chapter THREE Channel Reaction Rates 59 3.1 Overview In Chapter 2 we examined and discussed several approximations that could be used to make a model of CICR simpler and computationally less expensive. We dis- cussed the implications of these approximations and justified their use by analyzing their accuracy in modeling CICR compared to more traditional, more exact choices. Conversely, in this chapter, we aim to make CICR models more precise by examining closely one key aspect of the model: the channel reaction rates. One of the most critical choices to make in modeling calcium channels is the set of binding and dissociation rates of molecules to the receptor sites at a channel. These rates have a dramatic impact on the channel dynamics since the ratios of these quantities determine the relative amount of time spent in each state as well as, perhaps most importantly, the driving factors for the opening or closing of a channel. Furthermore, there appears to be very little agreement in the literature about the values of these transition rates (Table 3.2), which reflects little collective understand- ing of the interactions of individual molecules and ions with the channel receptors. For most of the transition rates, the literature does not agree even on the order of magnitude of these values. We therefore attempted, following the work of Swillens et al. (1999), to derive values for these rates based on experimental evidence, adapt- ing their method for a tetrameric channel structure. In addition to providing a starting point for constructing models of CICR, this method for selecting the rate constants may be used to probe more precisely the inner workings of the calcium channel receptor reaction dynamics. 60 3.2 Illustration In this section, we will examine the choices made for the values of binding and dissociation rates in one particular model, as an illustration of the importance of selecting these values carefully. Skupin et al. (2010) used the following transition rates (table recreated from sup- porting text for (Skupin et al., 2010)): a1 20 µM−1 s−1 IP3 binding with no inhibiting Ca2+ bound b1 20 s−1 IP3 dissociation with no inhibiting Ca2+ bound a2 .001 µM−1 s−1 Ca2+ binding to inhibiting site with IP3 bound b2 .03 s−1 Ca2+ dissociation from inh. site with IP3 bound a3 2.6 µM−1 s−1 IP3 binding with inhibiting Ca2+ bound b3 20 s−1 IP3 dissociation with inhibiting Ca2+ bound a4 .025 µM−1 s−1 Ca2+ binding to inhibiting site with no IP3 bound b4 .1 s−1 Ca2+ dissociation from inh. site with no IP3 bound a5 10 µM−1 s−1 Ca2+ binding to activating site b5 1.225 s−1 Ca2+ dissociation from activating site Table 3.1: Channel transition rates from Skupin et al. (2010). Refer to Figure 3.1 for a schematic illustrating the behavior of channel transitions. This model used an IP3 concentration of .1 µM and a base calcium concentration of .05 µM. There are several important implications to the stochastic channel dynam- ics. The most notable, perhaps, is that with an IP3 concentration of [IP3 ] = .1 µM, the binding rate for an uninhibited subunit is only 2 s−1 compared to a 20 s−1 disso- ciation rate. Consequently, the subunit has only a 9% chance of having IP3 bound at any given moment, and this is independent of any of the channel or cell dynam- ics. Because a tetrameric channel requires the concerted activation of three or more subunits in order to open, the probability of the channel having sufficient IP3 bound to be capable of opening is then only 0.28%, or about 1 in 350. There are gener- ally around 10-30 channels in each cluster, which means that having more than one channel in a cluster simultaneously open is rare with this model. 61 Additionally, even when all the chan- nels are closed and the cell is at equi- librium, the binding rate of calcium to the activating site is .5 s−1 compared to a dissociation rate of 1.2 s−1 , so when the cell is at equilibrium, there’s a roughly 40% chance of a channel subunit having activating calcium bound - four times more likely than the same proba- bility for IP3 . When the channel opens, the probability of having activating cal- cium bound rises to 90+%, but the IP3 distribution remains the same, which means that a channel opening doesn’t Figure 3.1: The original De Young-Keizer dramatically influence the firing prob- model. The quantities aj are binding rates ability of itself or other nearby chan- and the bj are dissociation rates. Fig- ure recreated from (De Young and Keizer, nels; this model has extremely weak in- 1992). (Repeated from Chapter 1.) terchannel coupling. Because the influence of a single channel decays rapidly upon closure of the chan- nel, if the channels themselves have weak coupling, then puffs (and therefore spikes) can only occur by coincidental simultaneous opening of multiple channels, which runs counter to the common intuition that a channel opening should provoke other channels to open. 62 (a) (b) (c) (d) (e) a1 20 400 80 20 ki+ 35.8 b1 20 52 .64 52 ki− 5 a2 .001 .2 .04 .001 kd+ .088 b2 .03 .2098 .48 3.077·10−4 kd− 1.67 a3 2.6 400 80 20 ki+ 35.8 b3 20 377 40 52 ki− 5 a4 .025 .2 .4 .025 kd+ .088 b4 .1 .0289 .0768 .0025 kd− 1.67 a5 10 20 15 10 ka+ 166 b5 1.225 1.647 12 122.5 ka− 65.9 Table 3.2: Comparison of channel rate constants from various models: (a) Skupin et al. (2010), (b) De Young and Keizer (1992), (c) R¨ udiger et al. (2007), (d) Thurley and Falcke (2011), (e) those computed in this chapter. As before (Tables 1.1, 3.1), the ai are binding rates and have units µM−1 s−1 and the bi are dissociation rates with units s−1 . Note: many articles list values of ai and di where di ≡ abii . For these articles, we have instead computed and listed the values of bi for easier comparison. We have included the corresponding nomenclature used in this chapter for convenience. 3.3 Construction Introduction In this section, we will look at data from experiments on calcium channels, and use this data to derive values for the channel transition constants. For this construction, we assume the usual eight-state DeYoung-Keizer model, and we assume precisely zero cooperativity. That is, one receptor being bound neither directly helps nor directly inhibits the binding of other receptors. In the context of the above table, this would mean that a1 = a3 , b1 = b3 , a2 = a4 , and b2 = b4 . Therefore, we seek values for six transition rates in total, which describe mutually independent stochastic events. Following the nomenclature of Swillens et al. (1999), we label these rates ki± for binding and dissociation of IP3 , ka± for binding and dissociation of Ca2+ at the activating site, and kd± for binding and dissociation of Ca2+ at the 63 deactivating (inhibiting) site. Table 3.2(e) lists the values for these rates that we will construct in this section, with the nomenclature used in this chapter included for convenience. Probability of opening We start with the assumption that ki− = 5 s−1 (Hannaert-Merah et al., 1995), which we note is well within the range of values for b1 in Table 3.2. We can then compute ki+ from a consideration of puff amplitude: if we take the number of chan- nels in a cluster to be 25, we take the size of a puff to be 5 times that of a blip, and we assume that most cluster events are puffs (that is, we assume strong coupling between channels in a cluster), then it follows that at any given moment, we expect on average one fifth of the channels to have IP3 bound to three or more subunits (so that the channel will open if flooded with calcium), and so each channel should have 1 a probability of 5 of being in such a state. The models Swillens uses does not have a subunit structure (the channel is monomeric), so the computation of ki+ is straightforward: the aforementioned prob- −1 ki+ [IP3 ] ki−  Ki ability is precisely k+ [IP ]+k − = 1 + [IP3 ] , where K i ≡ k+ is the dissociation i 3 i i 1 constant for this reaction. We can then set this quantity equal to 5 and easily solve for ki+ . With a tetrameric structure, the computation is only slightly more complicated. We define the following probabilities: pa is the probability that the channel has sufficiently many IP3 -bound receptors to be open (contingent upon cal- cium binding), pn is the probability that precisely n subunits are IP3 -bound, and  −1 Ki p = 1 + [IP3 ] is the probability for one specific subunit being bound. Then we have:   4 n pn = p (1 − p)4−n n 64 and so pa = p 3 + p4 = 4p3 (1 − p) + p4 = 4p3 − 3p4 Taking pa = 51 , we find p = .41755 (of the remaining roots, two are complex and the other is greater than 1). Setting  −1 Ki 1+ = .41755 [IP3 ] we then determine a value of ki+ = 35.8 µM−1 s−1 for a tetrameric channel with a cellular IP3 concentration of .1 µM. Mean open time We can determine two more rate constants, namely ka− and kd+ , from consideration of the mean open time of a channel. There are three ways an active subunit can deactivate, causing the channel to close: dissociation of IP3 , dissociation of Ca2+ from the activating receptor, and binding of Ca2+ to the inhibiting receptor. The expected rate for a subunit deactivation is the sum of the rates of these three processes, i.e., ki− +ka− +kd+ [Ca2+ ], and so the rate for channel closing is then 3(ki− +ka− +kd+ [Ca2+ ]), assuming that precisely three subunits are active and that a fourth subunit does not activate before a currently active subunit deactivates.1 Consequently, the mean open 1 Correcting for the possibility of four subunits being active would require a much more sophis- ticated approach that would necessarily involve all of the rate constants, due to the possibility of repeatedly bouncing between three and four active subunits. However, the probability of a channel having four simultaneously active subunits is bounded above by 3% (the probability of having IP3 bound to all four subunits), and so the correction to the rate constants should be small. 65 time for a channel can then be expressed as 1 − −1 Topen = ki + ka− + kd+ [Ca2+ ] 3 We already have a value for ki− (5 s−1 ), so two data points suffice to compute the other two rates. Swillens cites a patch clamp experiment (Bezprozvanny and Ehrlich, 1994) measuring mean channel open time as a function of cytosolic calcium concentration. With a concentration of [Ca2 +] = .2 µM, the mean open time was found to be 4.7 ms, and with a concentration of 500 µM, the mean open time was 2.9 ms. Using these values, we find ka− = 65.9 s−1 and kd+ = .088 µM−1 s−1 . Opening latency To determine ka+ , Swillens examined an experiment (Callamaras et al., 1998) wherein a cluster of calcium channels was suddenly flooded with a controlled con- centration of IP3 and the time needed for any channel to open was then measured. As the amount of IP3 introduced increased, the latency steadily decreased to a lim- iting value of about 80 ms. This implies that IP3 reaches saturation at a certain threshold, at which all channels have IP3 maximally bound, and channel opening is then gated by calcium binding to the activating receptors. For monomeric chan- nels, this once again leads to a very straightforward computation for ka+ . Assume a cluster with 25 channels and a base calcium level of .04 µM. Then, if IP3 is already bound, the expected time for the channel to open is the same as the expected time for activating calcium to bind: 1 tch = ka+ [Ca2+ ] 66 tch and the expected time for any channel in the cluster to fire is N where N is the number of channels. So in this case, ka+ = (25 · .08 s · .04 µM)−1 = 12.5 µM−1 s−1 . Adapting this for a tetrameric channel is dramatically more complicated both because the channel can take a lot of different paths through the state space to get to an open state, and because there is more than one state each channel can start in, as opposed to simply the “off” state in the above approach. To approach this problem, we define the channel as having five states corresponding to the number of active subunits (0 – 4), and we seek to compute the time of first arrival to the state 3 from an arbitrary starting state.2 For each state n, we define the following characteristic times: τn = time to transition away from state n to any other state τn0 = time to return to state n after transitioning to state n − 1 τn+ = time to first arrival at state n + 1 starting from state n Additionally, we define Pn to be the probability of the channel starting in state n, which depends only on the initial calcium concentration and the transition rates. With these definitions, the expected time to first channel opening is then τopen = P0 (τ0+ + τ1+ + τ2+ ) + P1 (τ1+ + τ2+ ) + P2 τ2+ First, τn is simply (sum of rates of transitions away from n)−1 , i.e., 1 τn = n· ka− + (4 − n) · ka+ [Ca2+ ] 2 We do not actually need to consider the state 4 because the “off” states - 0, 1, and 2 - cannot transition directly to 4 without going to 3 first. If the channel starts in state 3 or state 4, it will fire immediately upon binding IP3 , and there will be no contribution to the average latency. 67 For brevity, for the rest of this section, we will take k + to mean ka+ [Ca2+ ] and k − to mean ka− , since the calcium concentration is constant for this problem and since we are only concerned with the dynamics of the activating calcium receptors. Then 1 τn = nk − + (4 − n)k + For the state 0, the only transition possible is to the state 1, so we have τ00 = 0 and τ0+ = τ0 . For the rest of the states, the time to return to the state n after transitioning backwards to n − 1 is the time to transition away from n plus the time of first arrival from n − 1 to n. That is, τn0 = τn + τn−1 + To compute τn+ , we need to consider that the channel may go back to n − 1 an arbitrarily large number of times before proceeding to n + 1. The probability of transitioning from n to n − 1 is nk − pn = − + = nk − τn nk + (4 − n)k and the probability of transitioning from n to n + 1 is 1 − pn = (4 − n)k + τn Suppose the channel goes back to n − 1 a total of j times before going to n + 1. The probability of this occurrence is (1 − pn )(pn )j 68 and the expected time for this process to elapse is τn + jτn0 Therefore, the expected first arrival time from n to n + 1 is given by ∞ X τn+ = (1 − pn )(pn )j (τn + jτn0 ) j=0 ∞ ∞ ! X X = (1 − pn ) τn (pn )j + τn0 j(pn )j j=0 j=0 pn τn0   τn = (1 − pn ) + 1 − pn (1 − pn )2 pn = τn + τ0 1 − pn n 1 pn + = τn + τn−1 1 − pn 1 − pn 1 nk − + = + + + τn−1 (4 − n)k (4 − n)k Thus, we have 1 τ0+ = 4k + k−   1 1 τ1+ = + 3k + 3k + 4k + k− k−    1 1 1 τ2+ = + + 2k + k + 3k + 3k + 4k + Analogously to before, when examining IP3 transitions, we have that Pn (the prob- ability at equilibrium of n subunits having calcium bound to the activating site) is given by n  4−n k+ k−   4 Pn = n k+ + k− k+ + k− 69 Thus, we have reduced our expression for τopen to one involving only k − and k + . Furthermore, we have already computed the value of k − . Recall: k − ≡ ka− = 65.9 s−1 The latency time for the experiment was reported as 80 ms, so if the cluster size is 25 channels, then the expected latency of a single channel should be τ = 25·80 ms = 2 s. Using these values for τ and k − , we can solve for k + readily, for instance by using a software algebra system. The resulting value for k + ≡ ka+ [Ca2+ ] was found to be 6.64 s−1 . Since the experiment in (Callamaras et al., 1998) took place with a base calcium concentration of .04 µM, it then follows that ka+ = 166 µM−1 s−1 . Recovery time It remains to find a value for kd− . The primary effect of this parameter is clear; binding of calcium to the inhibiting sites of channels in an active cluster provides a mechanism to turn off that cluster, so the rate of dissociation of calcium from these sites determines a “recovery time” - the amount of time after the cluster shuts down before it can become active again. Swillens argues that when a cluster ceases activity, most or all of the channels in the cluster are desensitized - that is, they have calcium bound to the inhibiting receptors. This has the important implication that, once a cluster has become de- sensitized, a channel should recover much more quickly than the expected activation time for a channel. If this were not the case, then individual channels could fire upon recovery while the rest of the cluster was still desensitized, and therefore there would be a greater proportion of blips than experimentally observed. If we assume that most subunits are desensitized when a puff ceases activity, 70 then a channel will resensitize when three subunits resensitize. Thus, the criterion of fast channel resensitization is that the time for three subunits to resensitize is much faster than the single-channel latency. Recall that the expected latency for a single channel is 2 s. Therefore, we have: 1 3kd−  = .5 s−1 2s so kd−  .167 s−1 . This gives us one bound on kd− . To find another, we simply note that a cluster ought to be able to become fully desensitized on a reasonable timescale. This means that at calcium concentrations typical for a puff, after sufficient time has elapsed for the states to be well mixed, we should find a high proportion of subunits in an inhibited state. In other words, at these concentration levels, the probability of having calcium bound to the inhibiting site of any subunit should be relatively high. Assuming a cluster consists of a 5 × 5 array of directly adjacent channels, the average pairwise distance between channels is roughly 5Rch = .03 µm. The value of the steady-state distribution at this distance from an open channel is 5.3 µM. If we take this as the average concentration on a channel due to an open channel anywhere in the cluster, then the average concentration due to a puff event (equivalent to five open channels) would be 26.5 µM, and the binding rate to an inhibiting site would then be kd+ [Ca2+ ] = (.088 µM−1 s−1 )(26.5 µM) = 2.33 s−1 . Taking kd− = 1.67 s−1 (i.e., an order of magnitude higher than the previously derived constraint), and assuming a well-mixed state during a puff event, we then find the probability of an individual subunit not having inhibiting calcium bound to be −1 [Ca2+ ]  p≡ 1+ = .426 Kd 71 and the probability of a channel being desensitized is then Prob(2 or more subunits inhibited) = 1 − Prob(0 or 1 subunits inhibited) = 1 − (4p3 − 3p4 ) = .790 4 so about 5 of channels are desensitized after a sufficient amount of time has elapsed during a puff, which seems reasonable. As a final check, the same probability during equilibrium when channels are closed (using a base calcium concentration of .05 µM) is only .0000413 - as expected, the probability of a channel being desensitized is vanishingly small during cell equi- librium. Therefore, we are justified in taking kd− = 1.67 s−1 . 72 Figure 3.2: Dynamics of a single cluster with transition rates for tetrameric channels derived from experimental data. 3.4 Results We have now succeeded in constructing a full set of stochastic channel transition rates (Table 3.2(e)) from consideration of experimental results concerning various significant expected times. To test the validity of these rates, we incorporated them into a model of a single cluster of channels. The model used is the same as that for Figure 2.10a except for the choice of channel reaction rates: the biological parameters are taken from (Swillens et al., 1999), and we model a single cluster in the center of a cell-sized domain with Dirichlet boundaries. Instantaneous equilibration is assumed, and the stochastic process is governed by the binding and dissociation rates constructed in this chapter. Figure 3.2 shows an excerpt of the evolution of the average cellular concentration 73 Figure 3.3: Dynamics of a single cluster with parameters from (Swillens et al., 1999). Repetition of Figure 2.10a. due to the activity of this cluster. Comparison with Figure 3.3 shows that these rates produce results consistent with those of other models. Chapter FOUR Conclusion 75 4.1 Summary We have in this work examined several approximations, simplifications, and meth- ods for modeling aspects of calcium-induced calcium release. To reiterate, our pri- mary results were: • The assumption that the calcium concentration immediately equilibrates to a steady-state spatial profile upon the opening or closing of a single channel is a useful and reasonable approximation. • The assumption that the concentrations of any buffer species do not depend on location or time is a useful and reasonable approximation. • It is possible to model CICR treating a cluster as a single, fundamental stochas- tic object and obtain qualitatively reasonable results. • The assumption that the single-channel flux is concentrated at a point is not a reasonable approximation and is of questionable usefulness. Instead, taking the flux to be uniform over the channel volume offers an approach with less arbitrariness that is neither computationally expensive nor difficult to imple- ment. • It is possible to construct channel receptor binding and dissociation rates for a cluster of channels using data from experimental results. The use of such a construction yields a working numerical model for a cluster of calcium channels. When combining these results together, we can construct a full-cell model that can be executed start-to-finish in a reasonable amount of time on a home computer. 76 4.2 Full model 4.2.1 Single channel The first step, after selecting values for the various biological parameters, is to construct the steady-state profile due to a single open channel. To accomplish this, we use a simple second-order finite difference scheme for the reaction-diffusion equa- tions for the calcium and buffer concentrations. Here, we use the assumption of constant buffers to reduce the system to a single equation. Furthermore, we assume that the cell boundaries are sufficiently far away that we can employ radial symmetry to reduce the number of computed dimensions to 1, and additionally, that we can use Dirichlet boundary conditions, c(Rcell ) = c0 . Rc We take the spatial step to be dr = 5 , where Rc is the channel radius, and as- sume a uniform flux within the channel radius. This manifests in the finite difference scheme as an extra term on the right-hand side for r ≤ Rc . We evolve the concentration in time using a fourth-order Runge-Kutta scheme. We evolve for a sufficiently long time - around 3 ms - for the spatial profile to reach equilibrium. We then store the resulting equilibrium profile to be recalled later. 4.2.2 Single cluster Once we have the single-channel equilibrium profile, the next step is to compute the statistics for a single cluster of channels. We take a cluster to be a square array of 25 channels, directly adjacent to one another. We keep track of the calcium con- centration at the center of each of the 25 channels. Each channel is tetrameric - that is, it is comprised of four identical subunits. 77 ki+ 35.8 µM−1 s−1 IP3 binding ki− 5 s−1 IP3 dissociation ka+ 166 µM−1 s−1 Ca2+ binding to activating site ka− 65.9 s−1 Ca2+ dissociation from activating site kd+ .088 µM−1 s−1 Ca2+ binding to inhibiting site kd− 1.67 s−1 Ca2+ dissociation from inhibiting site Table 4.1: Summary of channel transition rates derived in Chapter 3. Following the work of De Young and Keizer (1992), we consider each subunit of a channel to be in one of eight possible states, corresponding to the presence or absence of IP3 , activating Ca2+ , and inhibiting Ca2+ . The binding and dissociation of these species from the receptor sites is governed by the transition rates given in Table 4.1, which are those we derived in Chapter 3. We use Gillespie’s algorithm to determine the length of each timestep and com- pute the binding and dissociation reactions of the receptors throughout the cluster. After each event, we determine whether a channel has just opened or closed: if a channel has three or four subunits with IP3 bound and Ca2+ bound to the activating site but not to the inhibiting site, it is considered to be open. When a channel is open, we update the concentration at each channel using the distance between each pair of channels and the steady-state spatial profile we com- Rc puted for a single channel. If the distance is not an integer multiple of 5 , we linearly interpolate between the two closest data points for the steady-state profile. The con- centration at a channel is the sum of the effects of all currently open channels plus the base concentration level, as described in Section 2.2.4. We run this simulation for a sufficiently long time - several seconds - to get a large amount of data, and we then compute the relative frequencies and average lengths of time for having different numbers of channels simultaneously open as detailed in Section 2.5. We repeat this process for a range of values of c0 since, as demonstrated in Figure 4.1, the channel dynamics can vary widely as a function of c0 . We choose 78 Figure 4.1: The fraction of time a single channel is open as a function of local calcium concentration. Repetition of Figure 1.5. values of c0 on a logarithmic scale in order to have finer granularity at concentrations close to the base concentration, since the effect of one cluster on another is usually rather subtle. 4.2.3 The cell Once we have statistics for a cluster of channels, we can model a single cluster much more quickly (Fig. 2.10b) or model several clusters in a cell (Fig. 2.12). For the latter, we randomly distribute roughly 30-50 clusters uniformly in an area in the center of the cell, far enough from the cell boundaries that we can ignore their effect. The statistics of a cluster change only due to the base local calcium concentra- tion; the effect of an active cluster on itself is already accounted for in the cluster statistics, so we do not need to be careful about computing peak concentration at an active cluster as we did for an individual channel. Additionally, the equilibrium profile due to a source is, outside the flux region, independent of the flux distribution provided the total flux is the same (Fig. 4.2). For these reasons, we can treat the 79 (a) (b) Figure 4.2: A comparison of the spatial equilibrium profile for (a) a single channel with point source or uniform flux, (b) a cluster of channels with uniform flux over the volume of a channel vs. volume of a cluster. Repetition of Figures 2.8 and 2.11. clusters as points, and so we keep track of the local calcium concentration at the center of each cluster. We model the cell dynamics in the same way as we did the single-cluster dynam- ics: we use the cluster statistics together with the Gillespie algorithm to determine the length of each timestep and what happens at each timestep. The principle dif- ference is that when a cluster undergoes an event, it does not have an effect on itself (again, because this is already accounted for in the cluster statistics). We com- pute the effect on other clusters using the distance between them, the single-channel steady-state profile, and the number of open channels. Since we are treating the cluster as a point, the profile for a cluster with n open channels is c(r) − c0 ) + c0 c(r) = n(˜ where c˜(r) is the steady-state profile for a single channel. More explicitly, at each timestep, we use the Gillespie algorithm to compute the length of the timestep and the next event, which will be equivalent to a channel 80 opening or closing, and be reflected in the state of the affected cluster. The calcium concentrations at every other cluster is then updated to reflect the new number of open channels in the affected cluster. Tracking the number of open channels at each timestep along with total elapsed time allows us to examine the time evolution of the average calcium concentration in the cell. 81 4.3 Further work While we have succeeded in probing several key facets of modeling CICR, and have demonstrated the plausibility of a computationally inexpensive model for a cellular domain with many clusters, there is still much to be learned. Open questions include: • Is the continuum approximation a reasonable assumption? How does a contin- uous diffusion model of CICR compare to a molecular dynamics model? • Are there situations in which the assumptions of instantaneous equilibration and constant buffers break down? Are these situations biologically realistic? • The curves in Figure 2.8 (point source vs. uniform flux) don’t match as closely as might be desired, implying there may be significant error in the concen- tration at the channel opening. Why does this discrepancy arise? Can the implementation of the channel flux be improved further? • Is cooperativity important for channel binding and dissociation reactions? What role does it play in the overall picture? • Using clusters as the fundamental stochastic object instead of channels has so far resulted only in very rough approximations to results of other models. How can this approach be improved? Can it produce quantitatively useful results? Appendix A The Gillespie-R¨ udiger Algorithm 83 A.1 Motivation Individual chemical reactions taking place on a relatively small scale can pro- duce effects on a much larger scale. In calcium channel dynamics, a single molecule binding to or dissociating from a receptor can have a dramatic impact on the dif- fusion dynamics (by causing a channel to open or close). In order to adequately capture such a phenomenon, it is therefore desirable to model the individual reac- tions stochastically, and on a timescale where no more than one reaction occurs per stochastic timestep. On the other hand, the individual reactions can take place on a much shorter timescale than the large-scale effects they produce. Diffusion of calcium ions away from an open channel can be much slower than the random binding and dissociation of ions from channel receptors. When modeling reactions numerically, in order to re- duce computational cost, it is desirable to maximize the timestep used for modeling the stochastic process, especially when the stochastic process is coupled to a slower process in this way. A traditional Markov chain modeling a Poisson process does not achieve this since, by definition, most timesteps will not see anything happen. Gillespie (1977) devised an algorithm to model stochastic processes whereby each stochastic timestep sees precisely one reaction take place, striking a balance between these two aims. The algorithm does this by, instead of fixing a stochastic timestep and computing what happens during this timestep, computing each timestep prob- abilistically. Specifically, each timestep is chosen to be the time between the most recent reaction and the next one. The nature of the next reaction is then also deter- mined probabilistically. Mathematically, this is done by indexing the possible reactions i = 1, · · · , n and drawing two random numbers r1 , r2 ∈ U (0, 1). We seek the stochastic timestep τ 84 and the index of the next reaction j, which are determined by   1 1 τ = ln α r1 and j−1 j 1X 1X αi < r2 ≤ αi α i=1 α i=1 Pn where αi is the reaction rate for the ith reaction and α = i=1 αi . The computed values of τ and j are then used to evolve the system from time t0 to time t0 + τ : if the system is in state X at time t0 , and the reaction Rj takes the system from the state X to Rj (X), then the numbers τ and j tell us that the system remains in state X for t0 ≤ t < t0 + τ and is in state Rj (X) at time t = t0 + τ . The primary drawback to this algorithm is that it only considers reaction rates that are constant between reactions. In modeling calcium channel dynamics, diffu- sion dictates that the calcium concentration be a continuous function of time. Since molecular binding rates depend on the quantity of species present, this means that the rate of the binding of calcium to channel receptors must also depend continu- ously on time. Alfonsi et al. (2005), therefore, extended Gillespie’s algorithm to include reaction rates αi = αi (t) that depend on time. Alfonsi’s results for determining τ and j are as follows: ˆ t0 +τ   1 α(t)dt = ln t0 r1 and j−1 j 1 X 1 X αi (t0 + τ ) < r2 ≤ αi (t0 + τ ) α(t0 + τ ) i=1 α(t0 + τ ) i=1 where t0 is the time at which the most recent reaction occurred. R¨ udiger et al. (2007) 85 provide a method for numerically computing τ from this equation and adapt this approach specifically for the framework of calcium ion channels. R¨ udiger’s algorithm for computing τ is examined at the end of this appendix. In this appendix, we will start with a physical consideration of molecular reac- tions. We will then derive from first principles a probability distribution P (τ, j; t0 ) for the occurrence of the first reaction after time t0 . From P , we will then derive the results given above, and demonstrate that they permit an algorithm for computing τ and j in a way that can be implemented in numerical models of chemical reactions. We will first recreate the derivation of Gillespie (1977) and verify these results in the case that reaction rates are time-independent between reactions, and then show how to generalize to the time-dependent case, as is done by Alfonsi et al. (2005). 86 A.2 Setup Suppose we have a chemical system with n species of reactants, S1 , S2 , · · · , Sn , with respective numbers of molecules N1 , N2 , · · · , Nn , and that there are m different possible reactions, R1 , R2 , · · · , Rm , between these molecules, with the j th reaction given by Rj := aj,1 S1 + aj,2 S2 + · · · + aj,n Sn → products We wish to determine the reaction rates αj corresponding to each reaction. The probability that a reaction of type Rj will take place in the time interval (t, t + dt) is the probability that a particular set of Rj reactants react in that time interval times the number of possible different sets of Rj reactants. For example, suppose R1 is characterized by S1 + S2 → S3 The probability of one particular molecule of species S1 and one particular molecule of species S2 reacting in the time interval (t, t + dt) is given by c1 dt where c1 is a constant that depends only on the properties of the reactants and the temperature of the system. In a state of thermal equilibrium, c1 is truly constant. Details about the constant c1 may be found in (Gillespie, 1977). Now, the probability of any reaction of this type occurring is the reaction prob- ability of a single set of reactants multiplied by the number of possible reactant sets. The number of distinct reactant pairs in this case is simply N1 N2 , and so the probability of an R1 type reaction occurring in (t, t + dt) is given by N1 N2 c1 dt 87 Note that the product S3 has no bearing on the dynamics of the reaction. Consider now the reaction R2 characterized by 2S1 → S4 with reaction probability constant c2 . (As before, the product S4 is completely arbitrary.) The number of possible reactant pairs in this case is   N1 N1 (N1 − 1) = 2 2 and so the probability of an R2 type reaction occurring in the interval (t, t + dt) is given by 1 N1 (N1 − 1)c2 dt 2 In general, if the reaction Rj is characterized by aj,1 S1 + aj,2 S2 + · · · + aj,n Sn → products with reaction probability constant cj , the probability of this kind of reaction occur- ring in (t, t + dt) will be      N1 N2 Nn αj dt ≡ ··· cj dt aj,1 aj,2 aj,n Note that because the cj are constant and the Nj only change when there is a reac- tion (assuming a closed system), then the αj are constant between reactions. 88 A.3 The probability distribution We are interested in the probability distribution for the occurrence of the next reaction: the time at which it will take place, and the kind of reaction it will be. Hence we define1 P (τ, j; t) such that P (τ, j; t)dτ is the probability that the next reaction that occurs after time t will be in the interval (t + τ, t + τ + dτ ) and will be of type Rj , given species quantities N1 , N2 , · · · , Nn at time t. Note that if we assume only one reaction can occur in an infinitesimal time interval (t, t + dt), then P (0, j; t) = αj by definition. To compute P (τ, j; t), we first define P0 (τ ; t) to be the probability that no reaction occurs in (t, t + τ ). Then the probability that the next reaction after time t occurs in (t+τ, t+τ +dτ ), and is of type Rj is equal to the probability that no reaction occurs in (t, t + τ ) multiplied by the probability that an Rj reaction occurs in (t + τ, t + τ + dτ ), i.e., P (τ, j; t)dτ = P0 (τ ; t)P (0, j; t + τ )dτ = P0 (τ ; t)αj dτ Computing P0 (τ ; t) is straightforward. First we note that the probability of no reaction occurring in (t, t + dt) is m X P0 (t + dt; t) = 1 − αj dt ≡ 1 − αdt j=1 So then P0 (τ + dτ ; t) = P0 (τ ; t)P0 (τ + dτ ; τ ) = P0 (τ ; t) · (1 − αdτ ) 1 The parameter t is unimportant for now, since the rates αj are constant between reactions, but it will become useful later when we no longer take the αj to be constant and so we include it here for consistency. 89 Thus 1 P0 (τ + dτ ; t) − P0 (τ ; t) = −α P0 (τ ; t) dτ In the limit as dτ → 0, we therefore have d ln(P0 (τ ; t)) = −α dτ Integrating from 0 to τ , we have: P0 (τ ; t) = P0 (0; t)e−ατ But P0 (0; t) = 1 trivially, so we have P0 (τ ; t) = e−ατ And finally, P (τ, j; t) = αj e−ατ provided τ and j take on appropriate values. That is:   αj e−ατ for 0 ≤ τ < ∞, j = 1, 2, · · · , n  P (τ, j; t) =   0 otherwise and we have the probability density function for τ and j. 90 A.4 Sampling the variables: Gillespie’s algorithm Now that we have the probability density for τ and j, we wish to find a method for sampling these quantities. Achieving this is straightforward in this situation, since τ and j are independent with this distribution: we have P (τ, j; t) = f (τ )g(j) with   αe−ατ for 0 ≤ τ < ∞  f (τ ) =   0 otherwise  αj for j = 1, 2, · · · , n   α g(j) =  0  otherwise and so P is simply the joint distribution for independent exponential and discrete random variables. Sampling these kinds of distributions is an elementary result, but we reproduce the process here for the purpose of comparison with the later result involving time-dependence. For a one-dimensional, continuous distribution, such as f (τ ), we sample τ in the ´τ usual way: let F be the cumulative distribution function for τ : F (τ ) = −∞ f (s)ds. Then F takes on values from 0 to 1. Set r = F (τ ). Then τ = F −1 (r), with the existence of the inverse function ensured by the continuity and monotonicity of F . Now, simply sample r uniformly on (0, 1) and use this relation to get a value for τ . To verify that this gives the desired result, we compute the probability that a value of τ sampled in this manner lies in the interval (s, s + ds). By definition of the density function, we expect this probability to be f (s)ds. Now, the probability of τ lying in (s, s + ds) is, by construction, the same as the probability of r lying in (F (s), F (s + ds)). Since r is uniform on (0, 1), this is simply F (s + ds) − F (s) = F 0 (s)ds = f (s)ds 91 as desired. For the function, f (τ ) = αe−ατ , we have F (τ ) = 1 − e−ατ and 1 F −1 (r) = − ln(1 − r) α Since r is sampled uniformly on (0, 1), r is statistically equivalent to 1 − r, and so we can slightly simplify this result:   −1 1 1 F (r) = ln α r We follow a similar procedure for the discrete distribution, g(j). In the case of a discrete “density”, with g(j) = pj for j ∈ N and ∞ P i=0 pi = 1, the cumulative distribution function is given by j X G(j) = pi i=0 To sample this distribution, we take r uniformly on (0, 1) as before, and choose the value of j that satisfies G(j − 1) < r ≤ G(j) To verify that this gives the desired result, we compute the probability that a sampled value of j will equal i. We expect this probability to be pi . By definition, the probability that j = i is the same as the probability that G(j) = G(i), which by construction is the same as the probability that G(i − 1) < r ≤ G(i). Since r is 92 uniform on (0, 1), this is the same as i X i−1 X G(i) − G(i − 1) = pi0 − pi0 = pi i0 =0 i0 =0 as desired. αj Thus, for the function, g(j) = α , we sample r uniformly on (0, 1) and then pick j such that j−1 j 1X 1X αi < r ≤ αi α i=1 α i=1 To summarize thus far, we have computed the probability that the next reaction after time t will occur in the interval (t + τ, t + τ + dτ ) and will be of type Rj . This probability is given by P (τ, j; t)dτ where   αj e−ατ for 0 ≤ τ < ∞, j = 1, 2, · · · , n  P (τ, j; t) =   0 otherwise We have shown that to sample τ and j from this distribution, we can sample r1 and r2 uniformly on (0, 1) and then compute τ by   1 1 τ = ln α r1 and pick j such that j−i j 1X 1X αi < r2 ≤ αi α i=1 α i=1 This method is Gillespie’s algorithm. 93 A.5 Time dependence Thus far, everything has assumed a closed, well-mixed system wherein the num- bers N1 , N2 , · · · , Nn do not change between reactions and hence the αj can be con- sidered constant when computing details about the next reaction. Such a situation is not always possible, however. In particular, if there is a concentration gradient of any of the reactant species, an external source of reactant species, or any other method by which a species can enter or leave a volume of interest, the αj , which depend on the quantities of the reactants, can change continuously in time, even between reactions. In this section, we will assume time-dependent reaction rates, attempt to follow the same procedure as in the proceeding sections, and highlight the differences that arise and how they may be treated. A thorough treatment of the time-dependent cast may also be found in (Alfonsi et al., 2005). Consider the previous situation with the one alteration that the αj now depend on time, i.e., αj ≡ αj (t) and α ≡ α(t). Computing P (τ, j; t) is identical in most re- gards. However, in computing P0 (τ ; t), we have P (0, j; t) = αj (t), and so we obtain the equation d ln(P0 (τ ; t)) = −α(t + τ ) dτ which has the solution ´ t+τ P0 (τ ; t) = e− t α(s)ds and so then the density function P becomes  ´  αj (t + τ )e− tt+τ α(s)ds for 0 ≤ τ < ∞, j = 1, 2, · · · , n  P (τ, j; t) =   0 otherwise Obtaining this result required little deviation from the previous approach, but now we encounter a difficulty: τ and j are no longer independent random variables! 94 Before, we could just factor the exponential distribution for τ and the discrete dis- tribution for j and treat them separately, but now the discrete factor αj depends on τ as well. Physically, this is exactly what we expect; since the reaction rates are changing in time, so too should the relative probabilities of the different reaction types. We may not be able to separate P , but we can still do something similar: Sup- pose we have a joint density f (x, y) where the random variables X and Y are not ´∞ independent. We use the marginal density for X, fX (x) ≡ −∞ f (x, y)dy, and use the same inversion algorithm as before. That is, ˆ x r1 = FX (x) = fX (x0 )dx0 −∞ and so, picking r1 uniformly on (0, 1), we set x = FX−1 (r1 ) Next, we use the conditional density for y given x, f (x, y) fX=x (y) ≡ fX (x) and proceed similarly: ˆ y r2 = FX=x (y) = fX=x (y 0 )dy 0 −∞ and so, picking r2 uniformly on (0, 1), we set −1 y = FX=x (r2 ) 95 To verify that this gives the desired result, we compute the probability Prob(x0 < x < x0 + dx0 and y 0 < y < y 0 + dy 0 ) with x and y sampled as described. We expect this probability to be equal to f (x0 , y 0 )dx0 dy 0 . The above probability is of course the same as Prob(x0 < x < x0 + dx0 ) · Prob(y 0 < y < y 0 + dy 0 |x0 < x < x0 + dx0 ) which, by construction of r1 and r2 , is equal to Prob(FX (x0 ) < r1 < FX (x0 + dx0 )) · Prob(FX=x (y 0 ) < r2 < FX=x (y 0 + dy 0 )) Since r1 and r2 are sampled uniformly, this is equal to (FX (x0 + dx0 ) − FX (x0 )) · (FX=x (y 0 + dy 0 ) − FX=x (y 0 )) = fX (x0 )dx0 · fX=x (y 0 )dy 0 f (x0 , y 0 ) 0 0 = fX (x0 ) dx dy fX (x0 ) = f (x0 , y 0 )dx0 dy 0 as desired. If either or both of x and y take on discrete values, we can use the same approach, but with the discrete formulation as presented earlier. This will be illustrated in applying this technique to our particular distribution. Recall:  ´  αj (t + τ )e− tt+τ α(s)ds for 0 ≤ τ < ∞, j = 1, 2, · · · , n  P (τ, j; t) =   0 otherwise 96 We first sample τ and then, given the fixed value for τ , sample j. We have: m X fT (τ ) = P (τ, j; t) j=1 m X ´ t+τ = αj (t + τ )e− t α(s)ds j=1 ´ t+τ m X − α(s)ds = e t αj (t + τ ) j=1 ´ t+τ − α(s)ds = α(t + τ )e t so then ˆ τ FT (τ ) = fT (s)ds −∞ ´ t+τ = 1 − e− t α(s)ds and so, setting r1 = FT (τ ): ´ t+τ r1 = 1 − e− t α(s)ds ˆ t+τ α(s)ds = − ln(1 − r1 ) t or equivalently, ˆ t+τ   1 α(s)ds = ln t r1 In the case that α does not depend on time, the integral evaluates to ατ , and so this result reduces exactly to the time-independent result. We have derived the sampling equation for τ , but we have now encountered another difficulty: this equation is not easily solved for τ ; the integral in general will not be able to be evaluated analytically, and so the left-hand side of this equation will not be invertible. We need another method to find τ . 97 A.6 Solving the integral: R¨ udiger’s algorithm We have just shown that ˆ t+τ   1 α(s)ds = ln t r1 ´ t+τ and wish to solve this equation for τ . If we define g(τ ) = t α(s)ds and ξ = ln( r11 ), we want to solve g(τ ) = ξ for τ , but g is in general not a function which is easily invertible. However, g is necessarily continuous and nondecreasing (since α ≥ 0), so there exists2 a unique3 value of τ that solves this equation, and we can find it simply by increasing τ until g(τ ) = ξ. This can be done easily by converting the integral equation for g into the initial value problem g(τ ˙ ) = α(t + τ ); g(0) = 0 and numerically integrating forward until the desired value for τ is found. This method for computing τ is the crux of R¨ udiger’s adaptation of the Gillespie algo- rithm. Now that we have a method for computing τ , it remains to sample j conditionally 2 It is conceivable that g asymptotically approaches some value less than ξ, but such situations would require the reaction rates to decay to zero, and would therefore require at least one species of reactant to evacuate the domain of interest entirely. Such a situation would require a different modeling approach. 3 It is also conceivable that there may be be intervals on which g is constant, but on these intervals, the reaction rates again would necessarily all be zero. Even so, the number of such plateaus is necessarily countable, and so the set of them has measure zero. 98 given this value of τ . We have: P (τ, j; t) fT =τ (j) = fT (τ ) ´ t+τ αj (t + τ )e− t α(s)ds = ´ t+τ α(t + τ )e− t α(s)ds αj (t + τ ) = α(t + τ ) So, given a value for τ , sampling j is exactly the same as before: given r2 sampled uniformly from (0, 1), pick j such that j−1 j 1 X 1 X αi (t + τ ) < r2 ≤ αi (t + τ ) α(t + τ ) i=1 α(t + τ ) i=1 We have therefore demonstrated the Gillespie-R¨ udiger algorithm in both the time-independent case and the time-dependent case. Appendix B Analytic Steady-State Solution due to a Single Open Channel 100 B.1 Overview It is desirable to have an analytic solution for the steady-state profile of a single open channel as it offers both a method of computation that can be incorporated into different models, as well as a benchmark for comparison. In this appendix, we will offer a brief derivation for the steady-state profile due to a localized source, originally printed in (Neher, 1986). Assume the flux source has zero physical extent and is embedded in an infinite flat membrane. We define the channel’s location to be the origin, and allow calcium to diffuse through the entire half-space exterior to the membrane. Then the calcium concentration, as a function of distance from the channel, r, is given by q + − kD b r I e c(r) = c0 + 4πDF r in this half-space, where c0 is the base concentration of calcium when the channel is closed (and so the limiting concentration as r → ∞), I is the open channel current, F is the Faraday constant, and the remaining quantities are as previously defined. 101 B.2 Setup and reduction Recall the reaction-diffusion equation for c: ∂c X  X = D∇2 c + ki− (Bi − bi (r)) − ki+ bi (r)c(r) + φj (t)δ(r − rj ) ∂t i=s,m,e j For simplicity, for the rest of this appendix, we will assume only one species of buffer. Then: ∂c X = D∇2 c − k + b(r)c(r) + k − (B − b(r)) + φj (t)δ(r − rj ) ∂t j When all channels are closed, there exists an equilibrium concentration c0 given by k + b0 c0 = k − (B − b0 ) or B − b0 c0 = K b0 k− where b0 is the buffer concentration at equilibrium and K ≡ k+ is the dissociation constant for the calcium-buffer reaction. (In actuality, models often specify B and c0 first, since those are the variables specified in experiments, and then let b0 be a derived quantity, but this is unimportant here.) Since we desire a steady-state profile for a single channel, we naturally seek a full solution of the form c(r, t) = c(r) = c0 + f (r). The angular independence arises from a simple consideration of geometry: the two simplest models involve either a channel far from any boundaries (in which case full radial symmetry is immediately apparent) or a channel embedded in a membrane. The membrane can be locally approximated as flat, in which case symmetry dictates that locally, the solution will be everywhere twice the “free channel” solution, and thus still radially symmetric in the half-space on the outside of the membrane. Substituting this ansatz into the full 102 equation yields a reduced equation for f : d2 f   2 df D 2 + = k + b0 f dr r dr d2 f 2 df k + b0 2 + = f ≡ µ2 f dr r dr D q k + b0 This equation has the two solutions 1r e−µr and 1r eµr . Since µ = D > 0, the latter is unphysical for large r and so we discard it. This leaves us with the full solution A −µr c(r) = c0 + e r where it remains to determine the value of A. 103 B.3 Flux condition We have determined that the steady-state solution has the form A −µr c(r) = c0 + e r q k + b0 where µ = D , and we seek to determine A. To do so, we consider a flux condition at r = 0. Suppose the open channel exhibits a steady current I. Then the molar flux per unit time out of the channel is given by φ= I 2F where F = 96485 C · mol−1 is the Faraday constant and the factor of 2 is due to the fact that a calcium ion has a charge of +2. Assume the geometry in which the channel is embedded in an infinite planar membrane, and consider a hemisphere of radius R centered on the channel. We integrate the steady-state equation over this space: ˚ ˚ ˚ 2 3 0= D∇ c dV + φδ (r) dV + (other terms) dV V V V ‹ ˚ 0=D ∇c · dS + φ + (other terms) dV ∂V V ˚ 2 dc 0 = D · 2πr +φ+ (other terms) dV dr r=R V noting that there is no flux through the membrane. If we let R → 0, the “other terms” vanish, and we have 2 dc I D · 2πr = −φ = − dr r=0 2F or 2 dc I r =− dr r=0 4πDF 104 Applying this to our solution:  0 2 dc A −µr 2 r = r c0 + e dr r=0 r r=0   2 −µr A A = r (−e ) 2 + µ r r r=0 = −Ae−µr (1 + µr) r=0 = −A I Therefore A = 4πDF and we have q + − kD b r I e c(r) = c0 + 4πDF r and so we have successfully computed the steady-state solution. Bibliography A. Alfonsi, E. Cances, G. Turinici, B. D. Ventura, and W. Huisinga. Exact simulation of hybrid stochastic and deterministic models for biochemical systems. ESAIM Proc., 14:1–13, 2005. F. J. Amador, L. Kimlicka, P. B. Stathopulos, G. M. C. Gasmi-Seabrook, D. H. Maclennan, F. Van Petegem, and M. Ikura. Type 2 ryanodine receptor domain A contains a unique and dynamic -helix that transitions to a -strand in a mutant linked with a heritable cardiomyopathy. J Mol Biol, 425(21):4034–46, Nov 2013. doi: 10.1016/j.jmb.2013.08.015. I. Bezprozvanny and B. E. Ehrlich. Inositol (1,4,5)-trisphosphate (InsP3)-gated Ca channels from cerebellum: conduction properties for divalent cations and regula- tion by intraluminal calcium. J Gen Physiol, 104(5):821–56, Nov 1994. M. Bootman, E. Niggli, M. Berridge, and P. Lipp. Imaging the hierarchical Ca2+ signalling system in HeLa cells. J Physiol, 499 ( Pt 2):307–14, Mar 1997. N. Callamaras and I. Parker. Radial localization of inositol 1,4,5-trisphosphate- sensitive Ca2+ release sites in Xenopus oocytes resolved by axial confocal linescan imaging. J Gen Physiol, 113(2):199–213, Feb 1999. N. Callamaras, J. S. Marchant, X. P. Sun, and I. Parker. Activation and co-ordination of InsP3-mediated elementary Ca2+ events during global Ca2+ signals in Xenopus oocytes. J Physiol, 509 ( Pt 1):81–91, May 1998. Q. A. Dang and M. Ehrhardt. On Dirac delta sequences and their generating func- tions. Applied Mathematics Letters, 25(12):2385–2390, July 2012. G. W. De Young and J. Keizer. A single-pool inositol 1,4,5-trisphosphate-receptor- based model for agonist-stimulated oscillations in Ca2+ concentration. Proc Natl Acad Sci U S A, 89(20):9895–9, Oct 1992. G. Dupont, L. Combettes, G. S. Bird, and J. W. Putney. Calcium oscillations. Cold Spring Harb Perspect Biol, 3(3), Mar 2011. doi: 10.1101/cshperspect.a004226. D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions. J Phys Chem, 8:2340, 1977. Z. Hannaert-Merah, L. Combettes, J. F. Coquil, S. Swillens, J. P. Mauger, M. Claret, and P. Champeil. Characterization of the co-agonist effects of strontium and cal- cium on myo-inositol trisphosphate-dependent ion fluxes in cerebellar microsomes. Cell Calcium, 18(5):390–9, Nov 1995. 105 106 Y. Jiang, A. Lee, J. Chen, V. Ruta, M. Cadene, B. T. Chait, and R. MacKinnon. X- ray structure of a voltage-dependent K+ channel. Nature, 423(6935):33–41, May 2003. doi: 10.1038/nature01580. D. O. Mak, S. McBride, and J. K. Foskett. Inositol 1,4,5-tris-phosphate activation of inositol tris-phosphate receptor Ca2+ channel by ligand tuning of Ca2+ inhibition. Proc Natl Acad Sci U S A, 95(26):15821–5, Dec 1998. E. Neher. Concentration profiles of intracellular calcium in the presence of a diffusible chelator. In U. Heinemann, M. Klee, E. Neher, and W. Singer, editors, Calcium Electrogenesis and Neuronal Functioning, pages 80–96. Berlin: Springer, 1986. I. Parker and Y. Yao. Ca2+ transients associated with openings of inositol trisphosphate-gated channels in Xenopus oocytes. J Physiol, 491 ( Pt 3):663–8, Mar 1996. S. R¨ udiger, J. W. Shuai, W. Huisinga, C. Nagaiah, G. Warnecke, I. Parker, and M. Falcke. Hybrid stochastic and deterministic simulations of calcium blips. Bio- phys J, 93(6):1847–57, Sep 2007. doi: 10.1529/biophysj.106.099879. M. Sanson, C. Ingueneau, C. Vindis, J.-C. Thiers, Y. Glock, H. Rousseau, Y. Sawa, Y. Bando, Z. Mallat, R. Salvayre, and A. N`egre-Salvayre. Oxygen-regulated protein-150 prevents calcium homeostasis deregulation and apoptosis induced by oxidized LDL in vascular cells. Cell Death Differ, 15(8):1255–65, Aug 2008. doi: 10.1038/cdd.2008.36. A. Skupin, H. Kettenmann, and M. Falcke. Calcium signals driven by single channel noise. PLoS Comput Biol, 6(8), 2010. doi: 10.1371/journal.pcbi.1000870. G. D. Smith, J. Wagner, and J. Keizer. Validity of the rapid buffering approximation near a point source of calcium ions. Biophys J, 70(6):2527–39, Jun 1996. doi: 10.1016/S0006-3495(96)79824-7. X. P. Sun, N. Callamaras, J. S. Marchant, and I. Parker. A continuum of InsP3- mediated elementary Ca2+ signalling events in Xenopus oocytes. J Physiol, 509 ( Pt 1):67–80, May 1998. S. Swillens, P. Champeil, L. Combettes, and G. Dupont. Stochastic simulation of a single inositol 1,4,5-trisphosphate-sensitive Ca2+ channel reveals repetitive openings during ’blip-like’ Ca2+ transients. Cell Calcium, 23(5):291–302, May 1998. S. Swillens, G. Dupont, L. Combettes, and P. Champeil. From calcium blips to calcium puffs: theoretical analysis of the requirements for interchannel communi- cation. Proc Natl Acad Sci U S A, 96(24):13750–5, Nov 1999. D. Thomas, P. Lipp, M. J. Berridge, and M. D. Bootman. Hormone-evoked ele- mentary Ca2+ signals are not stereotypic, but reflect activation of different size channel clusters and variable recruitment of channels within a cluster. J Biol Chem, 273(42):27130–6, Oct 1998. 107 K. Thurley and M. Falcke. Derivation of Ca2+ signals from puff properties re- veals that pathway function is robust against cell variability but sensitive for con- trol. Proc Natl Acad Sci U S A, 108(1):427–32, Jan 2011. doi: 10.1073/pnas. 1008435108. H. M. Viola and L. C. Hool. Cross-talk between L-type Ca2+ channels and mitochondria. Clin Exp Pharmacol Physiol, 37(2):229–35, Feb 2010. doi: 10.1111/j.1440-1681.2009.05277.x. J. Watras, I. Bezprozvanny, and B. E. Ehrlich. Inositol 1,4,5-trisphosphate-gated channels in cerebellum: presence of multiple conductance states. J Neurosci, 11 (10):3239–45, Oct 1991.