Simulia Multiscale Tissue

Ulrich Dobramysl1 , Sten R¨Udiger2 , And Radek Erban1

Abstract. Intracellular calcium is regulated in part by the release of Ca2+ ions from the endo- plasmic reticulum via inositol-4,5-triphosphate receptor (IP3R) channels (among other possibilities such as RyR and L-type calcium channels). The resulting dynamics are highly diverse, lead to local calcium “puffs” as well as global waves propagating through cells, as observed in Xenopus oocytes, neurons, and other cell types. Local fluctuations in the number of calcium ions play a crucial role in the onset of these features. Previous modeling studies of calcium puffdynamics stemming from IP3R channels have predominantly focused on stochastic channel models coupled to deterministic diffusion of ions, thereby neglecting local fluctuations of the ion number. Tracking of individual ions is computationally difficult due to the scale separation in the Ca2+ concentration when channels are in the open or closed states. In this paper, a spatial multiscale model for investigating of the dynam- ics of puffs is presented. It couples Brownian motion (diffusion) of ions with a stochastic channel gating model. The model is used to analyze calcium puffstatistics. Concentration time traces as well as channel state information are studied. We identify the regime in which puffs can be found and develop a mean-field theory to extract the boundary of this regime. Puffs are only possible when the time scale of channel inhibition is sufficiently large. Implications for the understanding of puff generation and termination are discussed.

Key words. Intracellular calcium, calcium puffs, multiscale modeling, stochastic diffusion

Ams Subject Classifications. 65C35, 92C42

1. Introduction. Intracellular calcium plays a major role in many signaling pathways, and regulates enzymatic activity, gene expression , and neural activity [60, 36]. In order to control a variety of cell functions, cells control the local cyto- plasmic calcium concentration via the exchange of Ca2+ ions with the extracellular space and the release of Ca2+ ions through channels situated on the membrane of reservoirs, such as the endoplasmic and the sarcoplasmic reticulum. Modulation of the Ca2+ release regulates muscle contraction, pathway cross-talk and mitochondrial activity, and disruption of these processes is associated with various diseases such as early-onset Alzheimer’s [5, 57], heart failure [3, 2] as well as psychological conditions such as bipolar disorder and schizophrenia . Hence, detailed knowledge about the underlying processes governing this signaling mechanism is required to allow progress in our understanding of these diseases.

In this paper, we focus on the release of calcium ions from the endoplasmic retic- ulum (ER) via inositol-4,5-triphosphate receptor (IP3R) channels. Upon binding of Ca2+ to binding sites on its cytosolic part, a channel opens and calcium ions flow from the ER into the cytoplasm. Channels occur in clusters of 10 to 20 channels [12, 47].

The opening of a single channel usually triggers the release of Ca2+ from other chan- nels in the same cluster due to the increased Ca2+ concentration in their vicinity . ∗Submitted to Multiscale Modeling and Simulation (April 1, 2015). The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013)/ ERC grant agreement no 239870.

Leverhulme Trust for a Philip Leverhulme Prize.

Deutsche

Forschungsgemeinschaft (RU1660 and IRTG 1740).

Arxiv:1504.00146V3 [Q-Bio.Sc] 13 Apr 2016

This mechanism results in a highly localized increase of the cytosolic calcium concen- tration. These “puffs” of calcium ions have been detected and analyzed in experiments by using fluorescent calcium buffers [58, 23].

Traditional modeling approaches of intracellular calcium dynamics are based on deterministic macroscopic rate equations [27, 11], however the intrinsically random, erratic nature of calcium signals in many cells or cell domains necessitates an ap- proach going beyond the deterministic regime [59, 24, 31, 44].

Progress Has Been

made by recognizing the importance of number fluctuations in the binding to the channels [48, 18] and using hybrid models, where the deterministic calcium concen- tration is coupled to stochastic channel binding models . Recently, it was found that local fluctuations stemming from diffusive noise (i.e. noise originating from the random movement of ions) have a crucial influence in calcium dynamics for clusters of intracellular channels, particularly in the inter-puffwaiting time but also on single-channel equilibrium behavior . Stochastic effects were also investigated for L-type calcium channels and RyR channels in the dyadic cleft [34, 49, 28]. Other studies consider a Langevin equation governing the fraction of open channels in a cluster . However, tracking the exact diffusion of individual ions in the complete computational domain is computationally intensive.

Here, for the first time, we apply spatial stochastic multiscale methods to model the dynamics of calcium puffs including the release of Ca2+ from IP3R channels and track individual ion positions in order to accurately incorporate diffusive noise. We take into account both activating and inhibitory channel properties. We study the dynamics of the calcium concentration as a function of the ion binding affinities and explore the regime where there exist puffs, as well as the regime in which channels do not close after their first opening. The paper is organized as follows: In the next Section 2 we discuss the spatial stochastic model for diffusion and channel gating used throughout this study. Section 3 discusses the multiscale approach we employ in order to reduce the computational effort required to track single ions and hence make this study feasible. In Section 4 we present results on the statistics of puffs extracted from simulated time series data, and study the transition between perpetually-open channel clusters and the parameter regime in which puffs can be observed. In addition, we develop a mean-field model for channel dynamics and use it to extract the boundary between the regimes. Finally, we summarize our findings in Section 5.

2. Spatial stochastic model for intracellular Ca2+ release. The spatial extent of our computational model consists of the three-dimensional domain Ωwhich models a part of the intracellular space. Ca2+ ions are able to undergo free diffusion in Ω.

They bind to and dissociate from binding sites on the channels, which are positioned on a small area of the domain boundary, corresponding to the membrane of the ER. The domain geometry and boundary conditions are specified in Section 2.2.

In the following sub-sections we discuss the components of this model in detail. The parameter values used throughout this study can be found in Table 1. 2.1. Diffusion - Brownian dynamics. A versatile method for the simulation of particles in a solution is given by Brownian Dynamics (BD). Collisions of particles with solution molecules lead to overdamped dynamics and random forcing on a suffi- ciently long time scale . Assuming that there are Q(t) free ions in the simulation domain at time t, the equation for Brownian motion of ions is given by

Table 1

Parameter values for the spatial stochastic simulations, for physiologically relevant conditions . where D is the diffusion coefficient of ions in solution, Xj(t) ∈Ωdescribes the tra- jectory of the j-th ion, and Wj is a three-dimensional vector of independent Wiener processes.

This approach dramatically reduces the dimensionality of the problem compared to molecular dynamics approaches wherein the degrees of freedom of every participating molecule need to be taken into account. Nevertheless, the computa- tional load is still high compared to deterministic PDE-based approaches to diffusion.

There are a number of approaches for simulating (2.1) in the literature, ranging from discretization with a fixed time step to event-based methods [52, 37]. In this paper, we discretize time with the time step ∆t and use the Euler-Maruyama discretization of equation (2.1), i.e. the position of the j-th ion is updated according to

J = 1, 2, . . . , Q(T) ,

where ξj is a vector of three independent normally distributed random numbers with zero mean and unit variance. 2.2. Domain geometry and boundary conditions. In our simulations, the computational domain is given as cube Ω= [0, L]3 where the value of L is specified together with other parameters in Table 1. Ca2+ ions are able to undergo free diffusion in Ωwhich we simulate using (2.2). A cluster of nine IP3R channels is positioned in a 3 × 3 grid with grid constant ℓ= 0.15 µm, centered in the z = 0 plane, i.e. the

3

positions of nine channels in the cluster are given as



. We found no significant effects from varying the channel spacing ℓ, hence we chose this particular cluster configuration for ease of implementation. Ions bind to and dissociate from binding sites on the channels. The boundaries of the computational domain at x = 0, x = L, y = 0, y = L, and z = L are constant-concentration boundaries, hence they absorb and introduce ions, such that in the absence of open channels the concentration of Ca2+ in the computational domain is held at its equilibrium value, c0, on average. The boundary at z = 0 is reflective, corresponding to the membrane of the ER. Hence, ions can enter and leave the domain via the boundaries, in addition to being introduced through open channels.

If the channels are closed, then the average number of free ions in the computa- tional domain, ⟨Q(t)⟩, is equal to c0|Ω| where |Ω| = L3 is the volume of Ω. Using our parameter values (see Table 1), c0|Ω| ≈1.5 × 103. The average number of ions in the simulation domain, ⟨Q(t)⟩, increases when the channels are open (by around two orders of magnitude). In order to simulate the system over sufficiently long time intervals, we will use a multiscale approach, which we describe in detail in Section 3.

2.3. Stochastic channel binding model. The conformational changes be- tween the open and closed states of IP3R channels are controlled by the binding of Ca2+ to activating and inhibitory binding sites. These channels consist of 4 sub- units, that can be itself in an active or a neutral/inhibited state. Each subunit has three different binding sites: An activating binding site for Ca2+ ions, an inhibitory binding site for Ca2+ ions, and an IP3 binding site.

To accurately capture channel gating events, we employ a simplified DeYoung- Keizer model . Here, we disregard IP3 dynamics and consider the effects of IP3 only via their influence on the dissociation constant bi of Ca2+ ions from inhibitory binding sites. Free Ca2+ ions can bind to the activating and inhibitory sites, while bound ions can dissociate from occupied sites. Therefore, there are two reversible reactions in our model for each subunit of a channel:

(2.4)

The rates aa and ai describe the binding affinity of a Ca2+ ion to an activating and inhibitory binding site, respectively, while the off-rates ba and bi describe the corre- sponding dissociation reaction rates. The variables Sa and Si describe the binding site state and can take values of 0 and 1. Channels consist of four subunits, with each subunit having one activating and one inhibitory Ca2+-binding site; see Figure 1. A model subunit can then be in three distinct states: neutral (no binding site occupied), active (only the activating site occupied) and inhibited (if the inhibitory site is oc- cupied regardless of the state of the activating site). A channel then opens when at least three of its four subunits are in the active state.

2.4. Ion binding dynamics. In order to precisely capture channel opening and closing dynamics, we need to implement a reversible BD binding model of Ca2+ ions

Ba

Fig. 1. (a) Sketch of a channel, containing 4 subunits with an activating ion binding site (green square) and an inhibiting site (red circle). (b) State space of a single subunit. Ions can bind to the activating and the inhibiting sites. A subunit is active only if it has an occupied activating site as well as an empty inhibiting site.

to their corresponding binding sites. Ions become binding candidates as soon as they enter a half-sphere of radius ϱ around a binding site (the channels are positioned on the ER membrane, hence only the half-sphere above them is available for ion binding). They are then allowed to bind to the site with a probability Pλ per time step that they spend in the binding region [35, 16]. An ion bound to an activating site (resp. inhibitory site) is allowed to dissociate with a probability 1 −exp(−ba ∆t) (resp. 1 −exp(−bi ∆t)) and placed at a distance of σ (unbinding radius) from the binding site. We use the values of binding and unbinding radii, ϱ and σ, as given in Table 1. Their values are chosen to be reasonable, such that no overlaps occur between channels and that they are large enough such that the BD simulation time step can be chosen reasonably large. Then we can pre-calculate the remaining parameter Pλ (binding probability) before the start of simulations using the approach described in [35, Section 5].

2.5. Channel opening and calcium flux. Each channel has four subunits.

I

(t)), j = 1, 2, . . . , 9, k = 1, 2, 3, 4, be the state of the activating (resp. inhibiting) site of the k-th subunit of the j-th channel in the cluster. Then the number of active subunits of the j-th channel at time t is



.

If N J

a(t) ≥3, then the channel is considered to be in an open state. When a channel is open, new ions are introduced with a rate of 3.12 × 105 s−1 (corresponding to channel current IC = 0.1 pA [8, 47, 51]) at the channel site, simulating the flux of Ca2+ ions out of an active channel.

The positions of released ions then evolve according to equation (2.2). 3. Multiscale approach. During a puffevent, a large number of ions are re- leased into the cytoplasm. If all channels are open continuously, then one can estimate, for the parameter values given in Table 1, that our computational domain may con- tain of the order of 105 ions during the peak of a puff. Since we are interested in time scales of minutes, tracking the individual position of this number of ions via Brownian dynamics becomes infeasible. However, high accuracy and individual ion

Ωbd

Fig. 2. Sketch of the computational domain Ω. It consists of two regions: The compartment- based part Ω\ ΩBD (everywhere except for the red box in the bottom center, with the compartment size h illustrated in the left bottom of the picture) and the BD domain ΩBD (shown as a red box in the bottom center). The boundary in both x- and y-directions as well as the boundary at z = L are constant-concentration boundaries, while the boundary at z = 0 is reflective. The channel cluster is positioned at the center of the z = 0 boundary (indicated as black squares), see (2.3).

positions are only needed in the vicinity of channel sites in order to ensure accurate implementation of BD binding dynamics, described in Section 2.4. Therefore, we split

× [0, Lbd]

containing the channel sites (red region in Figure 2), as well as the remaining space Ω\ ΩBD in which we will use a coarser description of ion movement as described in the next subsection. Here, LBD < L is the length of the edge of the cube ΩBD. We use LBD = L/5 in our simulations (see Table 1). In particular, we use BD simulations in a small fraction of 1/53 ≈0.8% of the computational domain Ω.

3.1. Compartment-based model for diffusion. We subdivide the region Ω\ ΩBD into compartments (cubes) with size h (illustrated in the bottom left of Figure 2) and employ a compartment-based method to simulate the movement of ions .

Ions are allowed to move between adjacent compartments with a rate of d = D/h2. The compartment-based algorithm only stores and evolves the number of ions in each compartment (rather than following individual ions). An event-based stochastic simulation algorithm is used to efficiently simulate this system. Several equivalent methods have been developed in the literature, such as the Gillespie algorithm , the Next Reaction Method , the Next Subvolume Method [14, 30], as well as the Optimized Direct Method .

In this paper, we employ the Next Reaction Method . For each possible move between two neighboring compartments, a putative time for the next jump of an ion

D Ac(T)

where t is the current simulation time, d = D/h2 is the jump rate per one ion, Ac(t) is the current number of ions in the compartment from which an ion is jumping and r is a random number uniformly distributed in the interval (0, 1). Clearly, the putative jump time (3.1) is infinity if Ac(t) = 0, i.e. if the corresponding compartment is empty. The putative times are smaller (on average) if the corresponding compartment contains more ions.

The putative jump times are inserted into a priority queue (a heap data struc- ture ), which enables us to efficiently extract the earliest jump time and thus the next jump. This move is then performed, and the numbers of ions in compart- ments (and the corresponding putative times and their entries in the priority queue) are updated. We then iterate this process by finding the minimal putative time and performing the corresponding ion jump at each iteration.

At the boundaries of the computational domain Ω(except for the boundary at z = 0, which is reflective), ions are absorbed (i.e. they leave the domain) or can enter with rates consistent with a constant equilibrium concentration of c0 outside the domain. The jump rate from outside the domain into a compartment just inside the domain boundary is c0Dh which is used, instead of dAC(t), in (3.1) to compute the corresponding putative times.

3.2. Coupling the BD simulation in ΩBD with the compartment-based approach in Ω\ΩBD. Several methods exist in order to couple the BD and compart- ment-based methods across their interface, such as the two-regime method (TRM) [19, 20] or the ghost-cell method . Here, we employ the TRM. In order to accurately capture diffusion across the interface, the jump rates from adjacent compartments into ΩBD need to be adjusted from the bulk rate d = D/h2 to the interface rate

D

π ∆t. These jump rates are used, instead of d, in (3.1) to compute the corresponding putative times. If the chosen move in the compartment regime is a jump from a compartment adjacent to the interface into ΩBD, the occupancy number of the compartment is reduced by one and a new ion is added in ΩBD at the distance x from the interface which is sampled from the probability distribution

,

where erfc is the complementary error function. The rate (3.2) and the distribution (3.3) are used to transfer ions from Ω\ ΩBD into ΩBD. In the opposite direction, the TRM transfers any ion from ΩBD which during the time step interacts with the interface. For a detailed discussion and the derivation of the probabilities and jump rates above, please see references [19, 20, 21].

We employ the particle-based simulation library package Tyche which implements the TRM . The TRM has also been recently implemented in Smoldyn .

Concentration [Μ M]

Fig. 3. Example concentration time trace for ai = 1 µm−1 s−1 and bi = 1 s−1. 4. Results. In the following section we present the results of our findings. Sim- ulations were performed with the parameter values listed in Table 1 unless noted differently. In Section 4.1, we discuss Ca2+ puffstatistics from simulation runs. We follow this by investigating the regimes in which puffs are visible in Section 4.2. A simple mean-field theory is then presented in Section 4.3. It describes the transition between puffs and perpetually open channels.

4.1. Puffstatistics. Figure 3 displays a sample time trace of the Ca2+ con- centration in the computational domain for ai = 1 µm−1 s−1 and bi = 1 s−1 for a time range of 20 s. The ion concentrations shows erratic calcium puffs with an am- plitude of up to 0.5 µm. The constant background of c0 = 0.02 µm corresponds to approximately 12 ions per µm3 and justifies our use of Brownian dynamics because fluctuations in such a small number of ions influence binding to the channels. The puffs are characterized by a sharp, almost instant increase in concentration, followed by an approximately exponential decay after the channel cluster closes. Puffs are sep- arated by a refractory phase in which channel subunits are inhibited and the cluster cannot open.

We analyze the concentration data given as the time values tj and concentration values cj, j = 1, 2, . . . , N, (where N is the index of the last data point of the simulation run) by choosing a puff-starting threshold concentration value

Var[Cj],

i.e. one standard deviation above the mean concentration ⟨cj⟩. For the data set shown below, the threshold was ct ≈0.1µM. The averaging in (4.1) is taken over all values of j, j = 1, 2, . . . , N, i.e. the sample mean and variance in (4.1) are estimated by

J=1

(cj −⟨cj⟩) . Puffs are identified by the concentration crossing the threshold point (4.1).

Puff

ending points are identified via a re-crossing of this value. The crossing points are

Given In The Ordered Index Set

S = {j | cj−1 < ct and cj ≥ct, j = 2, 3, . . . , N −1}

(Puffstarting Indices) And

E = {j | cj > ct and cj+1 ≤ct, j = 2, 3, . . . , N −1}

8

(puffending indices). The sets S and E are enumerated by the puffindex 1 ≤k ≤Np, where Np = |E| is the number of finished puffs in the data set, i.e. we disregard the last index in S if the last puffdid not finish and denote

E = {E1, E2, . . . , Enp},

where the elements in S (resp. E) are ordered, i.e. s1 < s2 < · · · < sNp (resp.

e1 < e2 < · · · < eNp). The inter-pufftimes are then given by

O ∩(Tp, ∞),

where the threshold tp = 0.25 s. This threshold filtering decreases the number of puffs considered by removing the puffs which are not well separated. Indeed, the inset in Figure 4(a) shows that the distribution of interpufftimes below the threshold follows an exponential distribution and is due to random channel reopenings stemming from residual calcium. The distribution of inter-pufftimes (histogram of set T ) is shown in Figure 4(a) for a set of 954 puffintervals. Its shape is similar to a Gamma distribution

1

Γ(κ)θκ τ κ−1e−τ/θ . Indeed, when choosing κ = 2.82 and θ = 0.64 s in equation (4.4) (such that the resulting distribution has the same mean and variance as our data), the comparison is good. As the Erlang distribution (the special case of a Gamma distribution with κ ∈N) results from summing κ exponentially-distributed random variables, we can speculate that an effective model of calcium channel opening might be described by a three-stage binding/unbinding process.

The Maximum Of The Inter-Pufftime

distribution, τm ≈1.81 s, can be interpreted as the typical time between puffs in this system. The distribution of puffamplitudes are determined by taking the maximum value



. The amplitude distribution (histogram of set A) is shown in Figure 4(b) and mirrors the broad distribution characterized in experiments .

The puffdurations are calculated by considering the full duration at half max- imum, which is a common criterion for determining the puffduration in the liter- ature . To this end, we calculate the indices of crossings of the half maximum

Rk

rk = max{j | cj−1 < ak/2 and cj ≥ak/2 and sk ≤j ≤ek},

Fk

fk = min{j | cj > ak/2 and cj+1 ≤ak/2 and sk ≤j ≤ek},

O

.

(D)

Fig. 4. Calcium puffstatistics for a set of Np = 1151 puffs from simulations with a total duration of 3 simulation hours. The mean and standard deviations of the distributions are indicated by the solid and dashed red lines, respectively. (a) Inter-pufftime distribution for 545 puffintervals after threshold filtering (4.3). The mean inter-pufftime is 1.81 ± 1.077 s. The green dashed line displays Gamma distribution (4.4) with κ = 2.82 and θ = 0.64 s. The inset shows the distribution of interpufftimes smaller than 250 ms with a mean of 0.59 ± 0.061 s. (b) Amplitude distribution.

The mean puffamplitude is 0.34 ± 0.12 µm (indicated by the red vertical lines). (c) Puffduration (full duration at half amplitude) distribution. The mean puffduration is 0.04 ± 0.02 s. (d) A two- dimensional histogram showing the strong correlation between the average number of open channels (over a puff’s duration) with the puffamplitude. The black line indicates a linear fit to the correlation data.

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.