Reactive clusters on a membrane

We investigate the reaction dynamics of diffusive molecules with immobile binding partners. The fixed reactants build clusters that are comprised of just a few tens of molecules, which leads to small cluster sizes. These molecules participate in the reaction only if they are activated. The dynamics of activation is mapped to a time-dependent size of an active region within the cluster. We focus on the deterministic description of the dynamics of a single cluster. The spatial setup accounts for one of the most important determinants of the dynamics of a cluster, i.e. diffusional transport of reaction partners towards or away from the active region of the cluster. We provide numerical and analytical evidence that diffusion influences decisively the dynamic regimes of the reactions. The application of our methods to intracellular Ca2+ dynamics shows that large local concentrations saturate the Ca2+ feedback to the channel state control. It eliminates oscillations depending on this feedback.


Introduction
A basic task of cells is to respond to inner and outer stimuli, which involves sequences of chemical reactions that form the signaling pathway. While numerous reactions take place between dissolved binding partners in the cytosol, others occur in the plasma or organelle membranes. A third class of reaction happens between cytosolic molecules and membrane-bound reaction partners. Often these reactions involve only a small region of the membrane at a time because the cell does not only interpret the total amount of additional substances, but also the precise location of its production. A well-known example of this type of reaction is in the formation of the second messenger cyclic adenosine monophosphate (cAMP). Dissolved adenosine triphosphate (ATP) binds to adenylyl cylcase, which is fixed in the plasma membrane, to synthesize cAMP. Depending on where this reaction occurs, it may trigger a break down of glycogen to glycose or gene expression [2].
The fixed membrane-bound reaction partners are often concentrated in a small membrane area, which we call a cluster. We will refer to the molecules in the cluster as fixed elements. Generally, not all molecules in the membrane patch participate in the reaction: only activated elements join in. Although the activating mechanism varies among membranebound processes, this principle is ubiquitous in biology. A single cluster with several tens of fixed elements is the focus of the present work.
The small number of fixed elements in a cluster implies small cluster areas. This strong localization of the reactions causes large spatial gradients of diffusing species. The cluster diameters are much smaller than the diffusion length of the species in solution. This imposes limits on the reactions due to diffusion of dissolved reaction partners towards or away from the cluster. Here, we will model the dynamics of an active cluster accounting explicitly for the diffusion of dissolved reaction partners.
Beside cAMP, cells employ other second messengers in the abovementioned signaling tasks. Ca 2+ is another prominent representative of a second messenger that participates in numerous processes [4]. It communicates the fertilization through an egg cell [20], controls apoptosis [25], and is vital for the excitation-contraction coupling in cardiac myocytes [5]. The mechanism by which a cell regulates the concentration of cytosolic Ca 2+ involves receptors on the membrane of intracellular compartments. Here, the inositol-1,4,5 trisphosphate (IP 3 ) receptor channel (IP 3 R) on the membrane of the endoplasmic reticulum (ER) serves as an illustrative example. Upon activation, the channel opens, which in turn results in a transient flux of Ca 2+ from the ER to the cytosol. Importantly, the open probability of the channel depends sensitively on the concentration of free cytosolic Ca 2+ . A moderate increase raises the tendency to release calcium, whereas a high concentration of Ca 2+ causes inhibition and closes the channel. Thus, the channel releases the species that controls its state.
The feedback of Ca 2+ on the channel dynamics becomes even more relevant when we take the spatial organization of IP 3 receptor channels into account. Generically, the receptor channels form clusters that are randomly distributed on the membrane of the ER. The typical inter-cluster distance is 2-7 µm [19]. The number of IP 3 Rs within a cluster has not yet been established experimentally. However, it is estimated that a cluster comprises 1-40 channels. Using the size of 18 nm for a single IP 3 R with all four subunits [12], we arrive at a cluster diameter of 18-100 nm. Thus, IP 3 receptor channels are tightly coupled by diffusion within a cluster because the Ca 2+ concentration decays on length scales of about 1 µm [26].
We see that intracellular Ca 2+ dynamics is another example of a reaction between partners fixed in a small membrane area and diffusing species. The channels are the fixed elements, Ca 2+ is the diffusing species and the reactions are production (here as release) and binding to and dissociation from the binding sites on the channel.
To date, modeling of intracellular Ca 2+ dynamics has proceeded along two distinct paths. In the first approach, deterministic models and spatially averaged concentrations have been used [7,15]. Since this approach neglects any spatial information such as gradients, the Ca 2+ concentration changes by only one order of magnitude during oscillations. However, simulations show that large gradients occur and that the Ca 2+ concentration changes by 3-4 orders of magnitude at an open cluster [26]. Hence, these models do not span the necessary range of Ca 2+ concentrations. In the second approach, stochastic methods have been applied taking the small number of fixed elements per cluster seriously [3,8,10,13,23]. For a review on intracellular calcium models, see [9]. The stochastic simulations have been performed with discrete sources, thereby incorporating gradients. They exhibit oscillations in regions of the parameter space where the deterministic behavior was non-oscillatory. To understand the loss of the oscillatory regime was one motivation of this study. However, our results will apply to many reactions with spatially localized reaction partners which create large gradients of diffusing species. Therefore we present a general formulation of the model in the following section and apply it to intracellular Ca 2+ dynamics in section 3.

General model
The aim of this section is to introduce a general formulation of a model that incorporates the areally restricted binding dynamics. As motivated in the introduction, we refered to the immobile molecules in the reacting region as fixed elements.
Each element occupies a non-negligible spatial fraction of a cluster because only several tens of closely packed elements reside therein. Therefore, the state of a cluster is well characterized by the area that is occupied by the activated elements. We call this region the active area of a cluster. It usually decomposes into unconnected parts. Nevertheless, we merge the area of all activated elements to one concentric patch. Its size equals the sum of the areas of all activated elements, and its radius is denoted by a. A change in the number of activated elements is translated into an alteration of the radius. If all elements are deactivated, a equals zero. a takes its maximum value a 0 when all elements are in the activated state. This mapping can be applied if the diffusion length of the diffusing species is larger than the cluster size. The procedure follows results by Swillens and Dupont [24].
We represent the cluster of fixed elements as a sphere with radius a. A typical value of the maximal radius a 0 is tens of nanometers. A second sphere surrounds the cluster concentrically. It has a radius b = 5-100 µm and describes the environment of the cluster.
The system contains m diffusive species. They are described by the concentration fields c(r, t) := {c 1 (r, t), . . . , c m (r, t)}. Their dynamics are of the general forṁ Here ∇ 2 r denotes the radial part of the Laplacian in three dimensions.
The functions f 1 and f 2 subsume the details of the dynamics for r < a and r > a, respectively. Most commonly, f 1 is dominated by production and f 2 by consumption.
The state of a fixed element is controlled by binding of the diffusive species to binding sites. The occupation of these binding sites determines the state. Usually, there are several binding sites per element. We denote with p i , i = 1, . . . , n the fraction of elements in the state i and refer to the p i as gating variables. The dynamics of the gating variables are governed by the general equatioṅ Some of the n states correspond to the activated state of the elements, so that the fraction of activated elements is determined by the sum of the corresponding gating variables. Consequently, the gating variables determine the radius a of the active area by some function f as a = f (c(a), p 1 (a), . . . , p n (a)).
We include a dependence on the diffusive species and on all n gating variables in f to account for the most general case. The values of the concentration fields and of the gating variables do not vary significantly within a cluster because the diffusion lengths are larger than a 0 . Therefore, we can pick a typical value to compute c and p i . We have chosen the value at the boundary of the cluster, which turns equation (3) into an implicit expression. The behavior of the equations (1)-(3) can be investigated by a bifurcation analysis, which determines the stationary states and their stability. We begin with the stationary states, i.e. the solutions of the equations The overbar indicates the stationary states of the Ca 2+ concentration profile and of the gating variables, respectively. The constantā denotes the stationary value of the active area. Equation (4a) can be treated separately for r <ā and r >ā due to the Heaviside step function. Since we demand the concentration profiles to be C 1 functions with respect to r, the matching conditions for the stationary solutions read as c i (ā) =c o (ā) and ∂c i /∂r(ā) = ∂c o /∂r(ā). The subscripts i and o indicate the inner and outer solution, respectively.
The computation ofā proceeds in two steps. Firstly, we solve equation (4a) with a fixed, but arbitrary value ofā. This results in a solution forc, which includesā as a still undetermined parameter. Secondly, we determineā self-consistently from equation (3) after inserting the solutions forc and forp i . Figure 2 shows a graphical method for determiningā. The dotted line is the bisection line, whereas the full lines represent the rhs of equation (3) for a specific model (see below). The stationary states are given by the intersections. Upon changing one parameter, the curve of f is shifted. It results in a change of the values or the number of fixed points. The existence of a saddle node bifurcation is easily deduced from such a plot. It occurs when f touches the bisection line. This is equivalent to the condition f (a) = 1.
Knowing the stationary points (c,ā), we investigate their stability. A linearization of the reaction-diffusion dynamics in (1) and the gating dynamics in (2) results in the equationṡ We define y as the perturbations δc of the diffusing species, and z i as the perturbations δp i of the gating variables. δ D denotes Dirac's delta function. Although a is not a dynamical variable in our model, it still changes in time. This is a consequence of equation (3) because a is computed from the evolving concentration fields and gating variables. a can be written as a =ā + δa with δa = δa(y, z). To evaluate δa from equation (3), we expand the expression a + δa = f (c(r) + y(r),p(r) + z(r))| r=ā+δa (6) to linear order: The derivatives of f have to be taken at (c(ā),p(ā)). The denominator only arises because of the evaluation of f at r = a in equation (3). When we combine y and z to an (n + m)-dimensional vector x = (y, z) t , the linearized equations take the matrix formẋ = Mx. If M can be diagonalized, the general solution for x is given by a linear combination of terms v i exp(ω i ). v i represents an eigenvector of M and ω i the corresponding eigenvalue. Consequently, the linear stability is uniquely determined by the eigenvalues of M. As shown in the appendix, M can be diagonalized. The eigenvalues λ i that originate from the gating variables constitute a subset {λ i } of all eigenvalues {ω i } and are all non-positive, i.e. λ i 0 ∀i, if the gating variable dynamics are rate equations derived from a master equation. Moreover, the eigenvectors v i of M that correspond to the eigenvalues λ i possess the structure v i = (0, . . . , 0, q i ) t with dim q i = n. This has two important consequences. Firstly, the eigenvalues from the gating dynamics do not contribute to any linear instability. Secondly, the solution for y does not depend on exp(λ i t). Therefore, the time dependence of y is solely governed by the eigenvalues that originate from equations (8) and (5a). Equation (5a) can be solved separately for r <ā and r >ā because of the Heaviside step function. The matching conditions are now y i (ā) = y o (ā), and the first derivative jumps according to due to Dirac's delta function in equation (5a). The continuity atā and equation (8)

Calcium dynamics
We now apply the method of section 2 to the dynamics of cytosolic calcium as a prototypic model system. The cytosolic Ca 2+ concentration c is governed bẏ The constants D and E denote the diffusion coefficient of Ca 2+ in the cytosol and the Ca 2+ concentration in the ER, respectively. The term k l (E − c) refers to a leak current, whereas k p c describes the calcium uptake by the ER through sarco-endoplasmic reticulum calcium ATPase (SERCA) pumps. Although it would be more realistic to model SERCAs using a Hill equation with coefficient 2, we approximate them by a linear expression for the sake of an analytical treatment. The last term in equation (9) describes the flux of Ca 2+ through IP 3 receptor channels. They represent the fixed elements, and the radius a of the active area is determined by the fraction of open IP 3 receptor channels. The state of a single IP 3 R is controlled by the state of its four subunits [6]. Each subunit expresses binding sites for Ca 2+ and IP 3 . Their occupation determines the state of the subunit. De Young and Keizer introduced a model to describe the dynamics of one subunit [7]. It consists of three binding sites: an activating Ca 2+ -binding site, an inhibitory Ca 2+ -binding site and an IP 3 binding site. Therefore, the state of a subunit can be specified by a binary triplet ij k. The first index represents the IP 3 binding site, the second the activating Ca 2+ -binding site and the last the inhibiting Ca 2+ -binding site. An index equals 1 when a site is occupied and 0 otherwise. The eight states that originate from the three binding sites are depicted in figure 1. Each state corresponds to one corner of the cube. The transition rates between the different states are indicated at the arrows. Binding of Ca 2+ and IP 3 is always proportional to the Ca 2+ and IP 3 concentration, respectively, whereas unbinding is independent of these concentrations.
Since a subunit is activated when IP 3 This is the specification of the function f (see equation (3)) for the Ca 2+ dynamics of the De Young-Keizer (DK) model. Our investigations of the Ca 2+ dynamics start with the stationary solution of (9) for an arbitrary value ofā, We applied ∂c/∂r(0) = 0 and c(b) = k l E/(k l + k p ) as boundary conditions. The latter complies with the base level of the system. The stationary value of p 110 in dependence on c and the IP 3 concentration I reads as Here, d 1 and d 3 denote the dissociation constants for IP 3 when Ca 2+ is not and is bound to the inhibiting site, respectively. The parameters d 2 and d 4 refer to the dissociation constants for the inhibiting Ca 2+ processes dependent on IP 3 binding. d 5 represents the dissociation constant for the activating Ca 2+ site [7]. The DK model assumes that the IP 3 dynamics is much faster than calcium binding and unbinding. This entails a fast equilibration between states with IP 3 bound and not bound. We eliminate the IP 3 dynamics adiabatically in the following and use the stationary values of the states with IP 3 bound and not bound. As shown in the appendix, the value ofp 110 is not changed by this approximation. Thus, the above analysis remains valid and we proceed to the stability of the fixed points.
The linearization of equation (9) results iṅ We define f c := k c (E−c)δa. Note that the inner concentration field y i is still restricted to r ā. In linear order, the varying value of a is translated into an additional flux density f c at the rim of the stationary active area. The solution of (14) is y(r, t) = exp(ωt)u(r) with and We used the boundary conditions ∂u/∂r(0) = 0 and u(b) = 0. The still unknown coefficients A and B are fixed by the continuity of u and the discontinuity of ∂u/∂r atā (see also equation (8)). The latter is a direct consequence of the last term in equation (14). The homogeneous system of equations for A and B possesses a non-trivial solution only if its determinant equals zero. This leads to the equation that determines the eigenvalue ω. η is given byāδa/sinh(k 1ā ), which can be cast into the form with and a 6 = (a 2 I + d 1 a 4 )/(I + d 1 ). If the system exhibits a zero eigenvalue bifurcation for a given pair (ā, I ), then ω = 0 should solve equation (17). Indeed, using the identity equation (17) can be transformed into 1 = f (ā). This is one of the conditions for a saddle node bifurcation.

Results and discussion
Some results obtained with the model concerning Ca 2+ dynamics have already been reported in [27]. Here, we consider parameter values beyond the range relevant for intracellular Ca 2+ release by IP 3 Rs to demonstrate model behavior which might be significant for other membranebound reactions. Diffusion of calcium plays a central role for the selection of dynamic regimes of the Ca 2+ dynamics besides the dynamics of the IP 3 receptor channel. Hence, we present results for D = 40 µm 2 s −1 and D = 220 µm 2 s −1 . Diffusion of Ca 2+ can be easily changed in experiments by application of Ca 2+binding proteins (buffers). The original DK model is based on a continuous distribution of IP 3 receptor channels. Two Hopf bifurcations bounding an oscillatory regime are the most prominent features of the Ca 2+ dynamics. We test whether this property is conserved when going from spatially continuous source terms in equation (9) to a discrete model. To this aim, we rescale the flux density with a typical cluster spacing R and a representative cluster radius a 0 while keeping the total flux constant, i.e. k c = k DK c R 3 a 3 0 . The resulting k c of 3 × 10 5 s −1 , which agrees well with realistic values [26], leads to a loss of the oscillatory regime. We find a single stationary state for all IP 3 concentrations, which is linearly stable. Decreasing the flux density by several orders of magnitude and thus approaching the Ca 2+ concentration values of the DK model does not restore oscillations. This holds because gradients still prevail.
These results do not mean that a model like equations (1)-(3) does not exhibit oscillations for some parameter values. In the following, we investigate oscillations of the model. In order to obtain oscillations, we choose parameter values supported by experiments, but a value of the activating Ca 2+ dissociation constant d 5 such that we obtain an oscillatory regime. We use d 2 = 3 µM for the inhibitory process in agreement with recent measurements [1] (Mak et al found dissociation constants up to 45 µM, but based on a specific model [18]). According to the experiments in [1], the coefficients for binding to the inhibiting site, a 2 and a 4 , are both set to 0.2 (µM s) −1 . The binding rate constant for the activating site can be evaluated from puff frequencies [28]. This implies a 5 1 (µM s) −1 . Assays of the dissociation constant for Ca 2+ activation yield values from 77 nM to 309 nM [16,17,22]. We chose d 5 = 0.823 µ M, which is motivated by the results depicted in figure 3. It shows the dynamic regimes of the model in dependence on d 5 and I. Oscillations occur only for larger values of d 5 .
Two saddle node bifurcations and a Hopf bifurcation terminate in a cusp. The oscillations arising at the Hopf bifurcation vanish via a bifurcation close to the lower saddle node bifurcation, which involves an increase in period. We assume it to be a homoclinic bifurcation. Typical oscillations are shown in figure 4 and in [27] for smaller values of a 5 . The pattern is the same for all examples. Upon increasing the IP 3 concentration, the system responds with a huge spike of Ca 2+ release and finally settles into small amplitude oscillations in the range of the inhibiting dissociation constants. The amplitude decays in space to negligible values within 1.6 µm (figure 4, right panel). Such oscillations had to be expected, since a feedback of Ca 2+ to a control process is only given in a range where the process is sensitive to changes of c. This range is around the dissociation constant d 5 for Ca 2+ activation and d 2 and d 4 for Ca 2+ inhibition. Changes of concentrations far above or below the dissociation constants do not exert a feedback on the dynamics [14]. Hence, the large difference between dissociation constants and concentration values occurring at the releasing channel leads to saturation of the control mechanisms of the IP 3 R by Ca 2+ . This is the reason for the small oscillatory regimes and amplitudes. Figure 5 shows that the structure of the bifurcation diagram does not change with a higher value of the activating dissociation constant. There is a single fixed point for almost all IP 3 concentrations. Stable limit cycles exist close to the bistable area. They only extend to IP 3 concentrations where the upper branch is unstable. Thus, we again find a very small band of IP 3 concentrations in which the system oscillates. The oscillations behave in the same way as described above.
The existence of oscillations does not solely depend on dissociation constants. It depends on rate constants as well. Therefore, we test the influence of the binding rate constant a for Ca 2+ activation and of the binding rate constants a 2 and a 4 for Ca 2+ inhibition on the stationary states. The concentration values of the stationary states are conserved, since we do not change the dissociation constants.
The Hopf bifurcation moves towards the left saddle node bifurcation when we increase the rates for the inhibitory processes. Figure 6 displays the difference between the IP 3 concentration values of the Hopf bifurcation and those of the left saddle node bifurcation. This difference decreases monotonically to zero with higher values of a 2 . Hence, the oscillatory regime shrinks for stronger inhibition. The opposite effect occurs for Ca 2+ activation. An increment of a 5 (while keeping d 5 constant) shifts the Hopf bifurcation to higher IP 3 concentrations.
Oscillations do not always disappear in a (putative) homoclinic bifurcation. The results in figure 7 illustrate a period doubling sequence while approaching the lower saddle node bifurcation. The left panel shows a period-2 example; the right panel depicts a period-4 example. Higher periods occur, too. Oscillations appear only when the upper branch is unstable for these parameter values, too, which leads to a small oscillatory regime. The oscillations are again considerably damped at a distance of 1.6 µm.
We find a different structure of the bifurcation diagram for parameter values like those in figure 8. For D = 40 µm 2 s −1 , there are two saddle node bifurcations and a Hopf bifurcation. This is similar to the results obtained above. However, increasing the diffusion coefficient changes the topology of

Conclusion and outlook
We have presented an extended study of a new modeling concept for diffusive species that react with immobile reaction partners. The fixed reactants are confined to small clusters. Our approach to describe the cluster dynamics is always applicable when the diffusion length is much larger than the cluster size. We applied the above method to the dynamics of intracellular calcium mediated by IP 3 receptor channels. The spatial restriction of the Ca 2+ flux to small membrane areas led to the disappearance of Ca 2+ oscillations computed in spatially continuous models. The enlarged values of the Ca 2+ concentration at the cluster resulted in a single linearly stable stationary state. Choosing smaller values of the channel flux constant k c did not restore Ca 2+ oscillations. Hence, the strong impact of spatial gradients on dynamic regimes will most likely also apply to localized reactions generating much smaller gradients than the gradients around a releasing Ca 2+ channel.
The flux constants we have used and which entail the large concentrations occurring at releasing clusters are based on recent simulation results. These simulations of Ca 2+ liberation close to experimental conditions show that Ca 2+ concentrations span the range of 25-170 µM in the center of an open cluster [26]. This is 3-4 orders of magnitude larger than the base level. At the same time, the concentration increases only 1-2 times base level at neighboring clusters. The existence of propagating waves proofs that the activation process is sensitive to these small concentration changes. Since the channel control processes experience concentration changes of several orders of magnitude and react to small changes already, the possibility of eliminating large concentrations from the dynamics by the rescaling of binding constants must be ruled out.
Even if the Ca 2+ concentration oscillated according to a deterministic model, these oscillations are not the oscillations observed in experiments [27]. Two observations led to this conclusion. Firstly, the range of IP 3 concentrations providing oscillations is too small. Secondly, the amplitude and the mean of the oscillations are already considerably damped in a distance of 1.6 µm from the cluster. Thus, they cannot represent the global Ca 2+ oscillations seen in experiments. The spatial damping of the oscillations applies to each reaction that produces a diffusing species.
Our findings support earlier results that deterministic models, including only activation by IP 3 , activation by Ca 2+ and inhibition by Ca 2+ , do not capture intracellular Ca 2+ oscillations [27]. Oscillations are driven by fluctuations in channel opening. The stochasticity of intracellular Ca 2+ dynamics is caused by the stochastic binding and unbinding of IP 3 and Ca 2+ to the small number of receptor molecules. Fluctuations cause spontaneous release in a single cluster. This leads to a release spike like the initial spikes in figures 4 and 7. Such a large amplitude event can lead to the opening of neighboring clusters and finally, via a nucleation process, to a wave traveling through the whole cell. If that occurs repeatedly, oscillation-like processes follow [8]. Thus, the amplitude of the initial spike is responsible for the amplitude of the oscillations. Nucleation may occur at different spots in the cell essentially at the same time when the IP 3 concentration is high [8], leading to almost regular periods.
Oscillations might as well be introduced by additional feedback, e.g., a Ca 2+ feedback on IP 3 production or the filling state of the endoplasmic reticulum. Our findings suggest that the initiation of global Ca 2+ release would still occur by wave nucleation, since the Ca 2+ dynamics would not undergo a local instability. The additional feedbacks would just modulate the nucleation probability periodically.
The present study sheds new light on the interplay between localization and fluctuations. This feature does not hold exclusively for IP 3 mediated Ca 2+ liberation, but other intracellular pathways should exhibit it, too. The involvement of a small number of proteins in cell signaling has already been established, see [11] for instance. It will be interesting to apply the current method to networks of strongly localized reaction sites. This will expand our knowledge on fluctuationinduced phenomena, which have been mostly studied by stochastic simulations. In addition, our results reveal that the deterministic limit still has surprises in store, as the period doubling (see figure 7) shows. Whether this sequence leads to chaos is the subject of ongoing investigations.
Moreover, we adopt a first-order scheme for the time integration and 50% of the stability criterion [21].

Glossary
Fixed element. Any reactive substance that is restricted to small spatial regions and that possesses two distinct states: activated and deactivated. A fixed element participates only in a reaction when it is activated.
Bistability. Existence of two linearly stable solutions of a nonlinear system.
Cytosol. The fluid portion of the cell in which all other internal compartments, e.g. the endoplasmic reticulum, are embedded.
Endoplasmic reticulum. An extensive membranous network within the cell, which serves as a major intracellular Ca 2+ store.
Hopf bifurcation. Emergence or disappearance of an oscillatory solution upon variation of a parameter in a nonlinear system.
Linear stability. Property of a stationary solution of a nonlinear system. The system relaxes back to this stationary state upon any infinitesimal perturbation from this state, if it is linearly stable, and amplifies any infinitesimal perturbation, if it is linearly unstable.

Receptor.
A specialized protein on a cell's membrane that binds to substances that effect the activities of the cell.
Saddle node bifurcation. Emergence or disappearance of two stationary states upon variation of a parameter in a nonlinear system. Second messenger. A chemical signal that relays a hormonal message from a cell's surface to its interior.