Effect of Buffers with Multiple Binding Sites on Calcium Waves

The existence and properties of intracellular waves of increased free cytoplasmic calcium concentration (calcium waves) are strongly affected by the binding and unbinding of calcium ions to a multitude of different buffers in the cell. These buffers can be mobile or immobile and, in general, have multiple binding sites that are not independent. Previous theoretical studies have focused on the case when each buffer molecule binds a single calcium ion. In this study, we analyze how calcium waves are affected by calcium buffers with two non-independent binding sites, and show that the interactions between the calcium binding sites can result in the emergence of new behaviors. In particular, for certain combinations of kinetic parameters, the profiles of buffer molecules with one calcium ion bound can be non-monotone.

a wide range of stimuli (typically mediated by the binding of hormones or neurotransmitters to receptors on the cell surface) initiate oscillations and waves of increased calcium concentration, and it is the dynamic properties of these oscillations and waves (in particular the period, localization, and amplitude) which act as the intracellular signal (Dupont et al. 2016;Falcke 2004Falcke -2005. Although the processes that control the cytoplasmic calcium concentration differ in detail between cell types, there is considerable overlap between the basic mechanisms, and thus it makes sense to talk of a generic calcium oscillation or wave model. Such a generic model relies on the excitable nature of the release of calcium from the endoplasmic reticulum (ER), whereby a small increase in cytoplasmic calcium concentration can lead to the release of much greater amounts of calcium from the ER, typically either through inositol trisphosphate receptors (e.g., in non-excitable cells), through ryanodine receptors (e.g., in striated muscle cells), or through both types of channels. Such excitable release of calcium from internal stores is generically called calcium-induced calcium release, or CICR (Dupont et al. 2016).
Because of the inherently excitable nature of CICR, an understanding of the dynamical behavior of calcium waves can be gained by the analytical study of the well-known FitzHugh-Nagumo (FHN) excitable model (Fitzhugh 1960(Fitzhugh , 1961Nagumo et al. 1962). Further, in many cell types calcium release happens on a much faster time scale than the recovery process, and thus the FHN model can be reduced to the single bistable model. This reduction is equivalent to the study of the leading front of wave solutions in the FHN model. However, the mechanism for controlling calcium waves differs from the traditional excitable mechanism in one crucial respect: the presence of large numbers of calcium buffer proteins. In normal conditions, more than 99% of cytoplasmic calcium ions are bound by buffer proteins. This is because free calcium ions are poisonous to the cell (because they activate so many things), and so buffers are used to control the concentration of free calcium. Therefore, in order to study calcium waves, it is necessary to investigate the dynamical behavior of buffered excitable systems. Along this line, numerous researchers have studied calcium wave propagation in buffered excitable systems with the presence of calcium buffers with a single binding site. The buffered bistable model was proposed by Sneyd et al. (1998). The existence and uniqueness of waves are established by a number of authors (Tsai and Sneyd 2005;Kaźmierczak and Volpert 2008a, b), and stability is shown in Tsai (2007). One needs to note that not every wave in mathematical sense is physiologically relevant. By this we mean that the elevated stable equilibrium of calcium concentration should expand into the region of lower stable equilibrium (ground state). A criterion for the existence of a physiologically relevant wave is derived in Tsai and Sneyd (2011, Proposition 3.2).
Previous work on the effects of calcium buffers on calcium dynamics has mostly assumed that calcium buffering can be modeled effectively by assuming that each buffer molecule has a single calcium binding site. This assumption is, in general, not satisfied; calcium buffers typically bind multiple calcium ions in a cooperative fashion, with the binding of one calcium ion affecting the rate at which calcium binds to the other binding sites (Schwaller 2010). However, recent work by Matveev and his colleagues (Matveev 2018;Chen and Matveev 2021) has extended the traditional analysis of buffers with a single binding site to include the effects of a second binding site, where the binding is cooperative. It is this general approach that we follow here.
It is important to note that cooperative binding is critical for our analysis. If the buffer molecules have several independent calcium binding sites we can represent the action of these buffers as a sum of buffer molecules each with a single binding site (see e.g., Sect. 6 in Kaźmierczak and Sneyd 2021) . However, the assumption of cooperative binding means that the binding sites are not independent, and the probability of calcium binding to any particular site will depend on which other binding sites are already occupied.
In particular, we consider the extreme situation in which the calcium ions can bind only consecutively, concentrating mainly on the case of two binding sites. Thus initially, calcium can bind only to the first site and the other sites are unavailable (hidden). After binding a calcium ion to this site, the second consecutive site is activated, e.g., exposed as a result of changing the buffer molecule conformation. This process can repeat consecutively. Unbinding of calcium follows an analogous process.
Another important fact is that due to the change of conformation, the buffer molecule can have a different diffusion coefficient than initially (Sorensen and Shea 1996), which additionally complicates the description. This phenomenon will be, however, not taken into account.
In this paper, we show that buffered systems with buffer molecules possessing multiple dependent calcium binding sites can have specific properties, not shared by the systems with buffer molecules having only one binding site (or independent binding sites).
The plan of this paper is as follows. In Sect. 2, we present the mathematical model. Section 3 is devoted to traveling wave problem and its fast buffering reduction. The validity of the fast buffering reduction is given in the appendix. Then, in Sect. 4, a criterion for the existence of a physiologically acceptable wave is derived. With this criterion, the effect of two-site buffers on the dynamical behavior of calcium waves is deduced in Sect. 5. Finally, conclusions and discussions are given in Sect. 6.

Buffered Bistable Model
Let us consider the case of buffer molecules which can bind m ≥ 2 calcium ions. For simplicity, we will confine to one single type of buffer molecules. Now set where M j represents the buffer molecule which binds exactly j calcium ions. Note that M 0 stands for the unbound form of buffers. Let the process of binding and unbinding of a calcium ion to the buffer molecule M j be described by the following reaction scheme with the kinetic constants k j + and k j − , respectively: The only solution of this equation which is bounded uniformly on IR is identically a constant. We thus have j=0,...,m where b 0 represents the constant which is the total concentration of buffer molecules at each spatial point. Consequently, we can replace M 0 by b 0 − j=1,...,m M j and consider the resulting system consisting of the first m + 1 equations of system (2). Then, system (2) is reduced to the following system It is seen that the equations for M j , j = 2, . . . , m, do not change, and that the equation for M 1 is modified. We note that for m ≥ 3 the obtained system is not monotone, so does not enjoy the comparison principle and the existence theorem of traveling waves in Volpert et al. (1994) due to the form of the equation for M 1 .

Buffered Bistable Model with Two Calcium-Binding Sites
To facilitate the discussion for the effect of multiple binding sites on traveling waves, we will only focus on the case where a single type of buffer is present, and such a type of buffer possesses exactly two calcium binding sites. Under the aforementioned assumption, system (4) is reduced to the following Buffered bistable system: A direct computation reveals that the constant states of system (5) are given by the following expressions: where the c j are defined in (3), and the M j 1 and M j 2 are given by

Traveling Wave Problem for Two-site Buffering Model
A traveling wave solution (c, M 1 , M 2 ) of system (5) connecting P 1 to P 3 is a solution of system (5) which is a function of the traveling wave coordinate variable ξ = x + vt, i.e.
and satisfies the boundary conditions where v is the wave speed.
We make one remark about the notation of the traveling wave coordinate. Since we use ξ = x + vt as the traveling wave coordinate, a wave solution with positive (resp. negative) wave speed v corresponds to a wave propagating from the right to the left (resp. from the left to the right).
Note that not all of wave solutions are biologically reasonable. When waves pass throughout the whole cytosol, the free cytosolic calcium is at a high stable steady state. Therefore, a positive wave speed v is required for a biologically reasonable wave solution.
In terms of the moving coordinate ξ = x + vt, the wave profile (c, M 1 , M 2 ) of a traveling wave solution of system (5) satisfies the following ordinary differential equations: subject to the boundary condition (9).

Rapid Buffering Model
In the following, we will apply the rapid buffering approximation (RBA) (Wagner and Keizer 1994;Keener and Sneyd 1998) to analyze traveling wave solutions to system (5). The RBA was first proposed by Wagner and Keizer (1994) to study the effect of rapid buffers with a single binding site on calcium diffusion and oscillations. Its central idea is to assume that buffering processes are very fast compared with the other reactions. Mathematically, this means that According to the RBA, we will use the second and the third equations in system (10) (the equation for M 1 ) to calculate the expression I b in the bracket of the first equation of system (10), which in turn reduce the full system to a single equation for the free calcium concentration. Here we give an outline of this approach. To proceed, the expression I b can be rewritten in the following form: where I b1 and I b2 denote the free terms in the second and third equation of system (10). The presence of factor 2 in the last equality of (11) follows from the fact that every M 2 molecule has two calcium ions bound. We will thus calculate I b1 and I b2 , which in turn gives the value of I b . Next, we rescale the expressions for I b1 and I b2 using the large parameter out of them. The rescaled expressions are then equated to zero, to obtain the dependence of the functions M 1 and M 2 on the variable c. Finally, by differentiation we calculate the D M M 1 − v M 1 and D M M 2 − v M 2 , and thus obtain the values of the unrescaled quantities I b1 and I b2 , and so the I b . Now, we carry out the RBA for system (10). The rigorous verification of this procedure is postponed to the appendix. As a large parameter L we take any of the quantities k 0,1 ± , e.g., L = k 0 + . Let us set Then, system (10) can be written as: As the coefficients κ are of the order of 1, so we can write I b1 and I b2 in the form: The quantities I 1 and I 2 do not depend on any large parameters, so, asymptotically, we are justified to demand I 1 = 0 and I 2 = 0. From the second equation we obtain Putting this relation into I 1 and equating it to zero, we obtain with K 0 , K 1 and K defined by Now differentiating (16) gives Hence, in the moving coordinate ξ = x + vt, we have Now, a further differentiation of M 1 (ξ ) with respect to ξ gives that Likewise, and with M 2 (ξ ) = θ 2 (c(ξ )) c (ξ ).
Finally, by plugging (11), (19), and (22) into the c-equation of system (13), it follows that the c-equation of system (13), and so that of system (10), can be rewritten as the following Rapid buffered bistable system: where It is worthwhile to note that Eq. (23) corresponds to the parabolic equation of the form: under the traveling wave ansatz c(x, t) = c(ξ ) = c(x + vt). Now, under the assumption of buffers with fast kinetics, a traveling wave solution of the Buffering system (5) can be approximated by a solution of the Rapid buffering system (23), as stated in the following proposition.
Proposition 1 (Rapid buffering reduction) Let scaling (12) be in force. Suppose that (v 0 , c 0 ) is a solution of the Rapid buffered bistable system (23) subject to the boundary conditions Then, there exists a large L > 0 such that for each L > L 0 , we can find a traveling wave solution (v L , c L , M 1,L , M 2,L ) of the Buffering bistable system (5) such that where M 1 (·) and M 2 (·) are defined by (16) and (20), respectively.
The validity of Proposition 1 will be shown in the Appendix. (See the proof of a more general Theorem A.1.)

A Useful Transformation
Let us note that given the function θ , the structure of Eq. (23) is the same as the structure of Eq. (12.34) in Keener and Sneyd (1998). Let us note that Eq. (23) can be written as where . (27) Now, as w (c) = D c + D M θ(c), then w is a monotonically increasing function of c. So, w and c are inverse functions of each other. Let us denote We have Equation (26) can thus be converted into an equation for w:

A Criterion for the Propagation of Buffered Waves
Theorem 1 (Existence and uniqueness of wave solutions of Rapid buffering system) Let w 1 and w 3 be the unique real numbers such that φ(w 1 ) = c b and φ(w 3 ) = 1 + c b .
Suppose that Then, there exist a positive increasing function c 0 : IR → IR and a positive number v 0 such that (v 0 , c 0 ) is a solution of the Rapid buffering system (23)-(24) subject to the boundary conditions Moreover, the solution pair The proof of the theorem uses the standard shooting method and the arguments from Sneyd et al. (1998) applied to the case of buffers with one binding site, but for completeness and better reference below, we insert it here.

Proof of Theorem 1
According to the analysis given in Sect. 3.4, it suffices to consider Eq. (31). To begin with, we will give local analysis of (31) around its singular points. To see this, we can write (31) as a first-order system: where The steady states of this system are equal to (w i , 0), i = 1, 2, 3, where Noting that φ and w are inverse functions of each other by (27) and (28), this implies Due to the fact that and using (30) and (34), we conclude that the states w 1 and w 3 are stable, whereas w 2 is unstable. In fact, we have The eigenvalues of the linearization matrices at the singular points (w i , 0), i = 1, 2, 3, are solutions to the equations and are given by the expressions: It follows that (w 2 , 0) is a repeller, whereas (w i , 0), i = 1, 2 are saddle points. The eigenvectors corresponding to the positive eigenvalue at (w 1 , 0) and to the negative eigenvalue at (w 3 , 0) can be written as: Before proceeding further, we make two observations. First, as d 1 < 0, d 3 < 0, then for v ≥ 0 the vector V 1 has a positive slope, while the vector V 2 has a negative slope. Moreover, given d 1 < 0, the function S 1 is a strictly increasing function of v with S 1 (d 1 , 0) = 2 √ |d 1 | and lim v→∞ S 1 (d 1 , v) = ∞. On the other hand, the function S 3 decreases with v ≥ 0. Next, if z(ξ ) = w (ξ ) > 0 for ξ in some interval (ξ 1 , ξ 2 ), then z can be treated also as a function of w via the identification z(w) = z(w(ξ )) = z(ξ ). Below, for simplicity, we will use the same symbol for z as a function of w. Now, we go to the shooting scheme. To proceed, let us note that, for given v ≥ 0, the trajectories of system (33 ) starting from (w 1 , 0) along the vector V 1 satisfy the equation Due to the form of f (·), we have f (φ(w)) < 0 for w ∈ (w 1 , w 2 ), and hence for v ≥ 0 , z(w) > 0 and z ,w (w) > 0 for w ∈ (w 1 , w 2 ]. Now, we consider the trajectory with v = 0. Indeed, from (38) the trajectory with v = 0 will lie above the w-axis for w ∈ (w 1 , w 2 ]. On the other hand, according to (37) and (32), for v = 0 the trajectory must touch the axis z = 0 for some w < w 3 . Taken together, it follows that for v = 0, the trajectory must intersect the axis z = 0 for some w 0 ∈ (w 2 , w 3 ). Next, we consider the trajectory with sufficiently large v. To proceed, let us first consider the case v > 0. Then, due to the fact that f (φ(w)) < 0 for w ∈ (w 1 , w 2 ), we have by (37): Since system (33) is autonomous, without losing generality, we can suppose that, given v > 0, w 2 = w(ξ )| ξ =0 . Thus, by the z-equation of (33), we conclude that, as long as It follows that, if v > 0 is sufficiently large, thenz(0) > z 0 /2, and so for all ξ > 0 As a result, given sufficiently large v > 0, the considered trajectory will cross the line w = w 3 for some z = z v > z 0 /2. The boundedness of z v is implied by the fact that for any solution z such that z(w) > 0 and z ,w (w) > 0 for all w ∈ (w 1 , w 3 ], we have where we used (32). This leads to the estimate To summarize, for v = 0 the trajectory of system (33) crosses the axis z = 0 for w ∈ (w 2 , w 3 ), and there exists v * > 0 sufficiently large such that the corresponding trajectory crosses the line w = w 3 for some finite z v * > 0. By using the continuity argument, we conclude that there must exist at least one v 0 ∈ (0, v * ), such that for v = v 0 the corresponding trajectory reaches the singular point (w 3 , 0). Moreover, such a v 0 is unique. Suppose to the contrary that there exist v 01 and v 02 > v 01 , for which the corresponding trajectories T 1 and T 2 join the points (w 1 , 0) and (w 3 , 0). Then, the trajectory T 2 starts at a bigger slope than T 1 , and, by the monotonicity of trajectories with respect to v ≥ 0, T 2 stays above T 1 for all w ∈ (w 1 , w 3 ). However, according to the first observation after (36) T 2 should cross T 1 from above, hence we arrive at contradiction. This completes the proof.

Separatrix Curve J C = 0 and Buffers' Effect
As it is seen from Theorem 1 and its proof, the 'sine qua non' condition for the existence of traveling waves with positive speed is the positivity of the integral However, as the value of the integral dc is known, similarly to Tsai (2013), we will consider only the integral We note that if J C (a, K 0 , K 1 ) > 0, then calcium traveling waves always propagate.
On the other hand, if J C (a, K 0 , K 1 ) < 0, then waves will cease to propagate if the product D M b 0 is large enough so that Therefore, the presence of buffers can prevent the propagation of calcium waves only if their kinetic characteristics satisfy J C (a, K 0 , K 1 ) < 0. This suggest a detailed study of the separatrix curve J C (a, K 0 , K 1 ) = 0, which will be done later. Before proceeding, let us note that, according to (1) and (17), the limit K 1 → ∞ corresponds to the case of buffers with single calcium binding site. Biologically, this can be explained by the fact that K 1 = k 1 − k 1 + , so K 1 tending to infinity means that the reaction Ca 2+ + M 1 runs only from the right to the left. Therefore, in fact there are no molecules with two calcium ions bound. Also note that in the limit K 1 → ∞, the right-hand side of (40) is reduced to Therefore, the condition that is a necessary condition enabling one-site buffers to slow down or stop the propagation of advancing calcium waves.

The Modified System Excitability Function a(K 0 , K 1 )
As we have noted, the separatrix curve J C (a, K 0 , K 1 ) = 0 plays a key role in the propagation of calcium traveling waves. In this subsection, we will investigate it in more detail. To proceed, we will use the equation J C (a, K 0 , K 1 ) = 0 to deduce that the excitability variable a can be represented as a function a(·, ·) of the kinetic characteristic pair (K 0 , K 1 ). Let us remind that according to (3) we have assumed the source function has the cubic form where a ∈ (0, 1) and S is a positive constant. Recall that the parameter a in the source function f characterizes the system excitability. Therefore, the function a(K 0 , K 1 ) characterizes the system excitability in the buffered system with the kinetic characteristic pair (K 0 , K 1 ). Further, for a given buffered system with the kinetic characteristic pair (K 0 , K 1 ), if 0 < a < a(K 0 , K 1 ), then J C (a, K 0 , K 1 ) > 0, and so waves always propagate. On the other hand, if a > a(K 0 , K 1 ), then J C (a, K 0 , K 1 ) < 0, and so waves may cease to propagate provided if the product (b 0 D M ) of the concentration of total buffers and their diffusion coefficients are large enough. These follow from the fact that the function J C (a, K 0 , K 1 ) is decreasing in a, as shown in the following lemma.
Proof The proof follows from the fact that due to (43) With the use of Lemma 1, the separatrix curve J C (a, K 0 , K 1 ) = 0 gives rise to the existence of the modified excitability function a(·, ·), as shown in the following lemma.
Lemma 2 (Existence of the modified system excitability function a(K 0 , K 1 )) There exists a function a : For all 1 > 0, 2 > 0, the function a(·, ·) is of C 1 class of their arguments on the set Proof The existence of the function a(·, ·) follows from the implicit function theorem and Lemma 1. Likewise, the differentiability follows from the fact that, for j = 0, 1, we have from which (and Theorem 1) follows the boundedness and continuity of the derivatives.
It should be noted that for It follows, according to (44), that a K 0 (0, K 1 ) → ∞ as K 1 → ∞. These properties of the function a(K 0 , K 1 ) are illustrated in the right panel of Fig. 2. Similar singularity analysis can be carried out for Though the exact expressions of the modified excitability function a(·, ·) can be obtained (e.g., by using Mathematica code), they very complicated, thus it is better to use numerical simulations for the analysis. To fix our attention, in the remainder of this section, we will choose for the calcium ground concentration, unless it is indicated differently. The resting cytoplasmic calcium concentration is typically around 100 nM, but can vary anywhere from 20 nM to over 200 nM (Clapham 2007). The value we choose here is thus slightly on the high side, but still within the physiological regime.
For the case of one-site buffers, the exact expression of the function a(·, ·) is given by Sneyd et al. (1998). For two-site buffers, approximate (asymptotic) formulae can be only derived in some specific cases. For example, this can be done for very large values of K 0 , as it is shown in Sect. 5.2.2.
The shapes of the curves J C (a, K 1 , K 1 ) = 0 for some values of a ∈ (0, 1) are depicted in Fig. 1. The critical point (K * 0 , K * 1 ) ≈ (0.133, 2.376) of the modified excitability function a(·, ·) for which ∂a/∂ K 0 (K * 0 , K * 1 ) = 0 and ∂a/∂ K 1 (K * 0 , K * 1 ) = 0 divides the family of the curves J C (a, K 1 , K 1 ) = 0 into two classes. Specifically, let a * = a(K * 0 , K * 1 ) ≈ 0.417425 be the critical excitability parameter. Then, for a ∈ (0, a * ), the curve J C (a, K 1 , K 1 ) = 0 consists of two parts: one is opening upwards, and the other is opening downwards. On the other hand, for a ∈ (a * , 1), one component of the curve J C (a, K 1 , K 1 ) = 0 opens to the left, and the other opens to the right. This can be seen from the left and right panels of Fig. 1. The presentation of the form of the curves J C (a, K 1 , K 1 ) = 0 is continued in Fig. 2. In the left panel of Fig. 2 we show the region of smaller values of K 1 , whereas in the middle panel the region with larger values of K 0 . Finally, in the right panel we show the form of the curves a(K 0 ; K 1 ) as a function of K 0 for some given values of K 1 .

The Effect of Two-site Buffers
In this subsection, we analyze of the effect of binding sites on traveling wave solutions.

Admissible Region A for Wave Propagation
In this subsection, we discuss the region in the K 0 K 1 plane for which the corresponding wave propagates from the right to the left, hence with v > 0. If these conditions are satisfied we will simply say, for brevity, that the waves propagate.
To proceed, recall that for a fixed a ∈ (0, 1), if the kinetic pair (K 0 , K 1 ) lies in the region then the corresponding wave always propagates. On the other hand, if (K 0 , K 1 ) falls outside the region A(K 0 , K 1 ), then waves cannot propagate provided if the product b 0 D M is large enough. In order to be compared with one-site buffers, we need to define a critical K 0 . Indeed, for a given a ∈ (0, 1/2), there exists a critical K s 0 ≥ 0 such that Therefore, for the case of one-site buffers, when K 0 > K s 0 , waves always propagate, whereas for K 0 ∈ (0, K s 0 ), waves will be stopped provided if the product b 0 D M is large enough. As it is seen from the left panel of Fig. 3, the minimal excitability parameter a for which the action of buffer molecules can stop the advancing wave propagation is approximately equal to 0.331. On the other hand, for two site buffers, it can be shown that the minimal value of a is approximately not bigger than 0.254. (For example, for K 0 = 10 and K 1 = 0.00007 it is equal to 0.254.) These differences are also expressed in the regions in the K 0 K 1 -space, in for which the advancing waves can propagate independently of the values of D M b 0 .
With the use of K s 0 , the region A is decomposed into two subregions A + and A − , defined by The regions A, A + , A − , for two characteristic values of the parameter a, are depicted in Fig. 3. As shown in Fig. 3, one of the boundaries of the region A − tends to the vertical We first discuss the case where (K 0 , K 1 ) ∈ A + . For one-site buffers, the corresponding wave always propagates. For two-site buffers, although the corresponding wave always propagates if K 1 is large enough. On the other hand, the wave can fail to propagate provided K 1 is small and the product b 0 D M is large enough. As shown in Fig. 3, such a region (bottom portion of purple color with K 0 ≥ K s 0 ) is relatively small. (Note that the graphs in the middle and right panel are in the logarithmic scale.) Hence, we can come to an approximate conclusion that two-site buffers retain propagation for the kinetic pair (K 0 , K 1 ) ∈ A + for which one-site buffered waves always propagate.
Next we consider the case where (K 0 , K 1 ) ∈ A − . For one-site buffers, the corresponding wave can be stopped provided b 0 D M is large enough. However, for two-site buffers, some of the corresponding waves can propagate. Therefore, we can conclude that two-site buffers can promote propagation for the kinetic pair (K 0 , K 1 ) ∈ A − for which one-site buffered waves can be stopped if the product b 0 D M is large enough.

Can Buffers with Multiple Binding Sites Facilitate Calcium Wave Propagation?
One can ask an alternative question: Can the presence of buffers facilitate the propagation of calcium waves? Let us consider the extreme case where waves do not exist in the system with the absence of buffers, that is, Then, our question is the following. Suppose that condition (46) holds. Can the addition of buffers promote the propagation of waves? Intuitively, this seems to be impossible. However, as we will see, one-site buffers cannot promote the propagation, whereas two-site buffers can empower the propagation of advancing waves.

Single binding site
For the case of one-site buffers, the answer is negative, as it follows from the form of the integral J C (a, K 0 ) given by (41) and representing (up to a positive factor) the influence of buffers. Indeed, using the mean-value theorem for integrals and the fact that f is one sign on the intervals (c b , c b + a) and (c b + a, c b + 1), we can find c − ∈ (c b , c b + a) and c + ∈ (c b + a, c b + 1) such that This in turn implies that Hence, in this case the waves cannot exist by Theorem 1.

Two Binding Sites
For the case of two-site buffers, a successful propagation depends on whether the quantity is positive or not, as indicated by Theorem 1. Since the integral 1+c b c b f (c)dc is negative, we need to look for a triplet (a, K 0 , K 1 ) such that J C (a, K 0 , K 1 ) > 0. Then, for such a triplet (a, K 0 , K 1 ), the corresponding wave will propagate provided the product D M b 0 is large enough.
First, the right-hand side of (40) can be written in the form where for f given by (43), It follows that for K 1 > 0, in the limit K 0 → ∞, For c b = 0.2 and K 1 → 0, a(∞, K 1 ) → 0.571. We have thus shown the following theorem.
For the excitability parameters a ∈ (0, 1/2) the statement of the theorem is obvious as the limit K 0 → ∞ corresponds to the situation where calcium ions are completely released from bound buffers by the reaction scheme (1) and the relation K 0 = k 0 − /k 0 + . On the other hand, for a > 1/2 the claim of the theorem seems paradoxical and it is hard to find its physical explanation. Theorem 2 is confirmed numerically. Their results are presented in Fig. 4. It is numerically shown that for K 1 = 1 (left panel) and K 1 = 0.1 (right panel), a(K 0 , K 1 ) > 0.53. For these parameters and a = 0.53, there exist traveling wave solutions with v > 0. Moreover, this result holds not only for traveling wave solutions to the asymptotic Eq. (25) (heteroclinic solutions to Eq. (23), but also to traveling waves of an initial system (5) representing Eq. (23) with sufficiently large L.

Monotonicity Versus Non-monotonicity of Wave Profiles
As the initial system of reaction diffusion Eq. (5) is in general a non-monotonic one, one can expect that the profiles of some of its components can be also non-monotone. The same remark refers to the profiles of c, M 1 and M 2 defined by system (23), (18) and (21). First, let us note that the profile of the c-component of the traveling wave solutions to Eq. (23) is always monotone increasing, as can be seen from Lemma 1 and its proof. More precisely, we have Therefore, it remains to consider the profiles M 1 and M 2 of bound buffers.

Single Binding Site
This case corresponds to K 1 = ∞, hence from (18) and the fact that the c-component of wave solutions is monotone increasing, it follows that the M 1component of wave solutions is monotone increasing. Therefore, for the case of one-site buffers, the profiles of the concentrations of free calcium ions and free/bound buffers of wave solutions are always monotonic.

Two Binding Sites
For the case of two-site buffers, the situation is more complicated. Let us consider the profile of the M 1 -component of wave solutions. For the M 1 -component of wave solutions to be increasing, by (49), we have that M 1,c (c(ξ )) = M 1,c c (ξ ) > 0 on IR. With the aid of (18), this is equivalent to the inequality K 0 − K −1 1 c 2 > 0 for c ∈ (c b , c b + 1), and hence that K 0 > K −1 1 (c b + 1) 2 . Therefore by (18), we can conclude that Similar arguments give that For each K 1 > 0 fixed, the above two observations motivate us to define the following two curves in the parameter (a, K 0 )-plane: Another feature of the profile of the M 1 -component of wave solutions is that in the presence of two-site buffers, it may not be heteroclinic. Indeed, for the profile of the M 1 component to be heteroclinic, by (16), we must have that M 1 (c(−∞)) = M 1 (c(∞)), and so that M 1 (c b ) = M 1 (c b + 1). This motivates us to define the curve m in the parameter (a, K 0 )-plane as follows: (the dashed (middle) horizontal line in Fig. 5 and 6) Therefore, for the parameter pair (a, K 0 ) lying on the curve m , the profile of the M 1 -component of the corresponding wave solution is homoclinic. Further, one can verify that for the parameter pair (a, K 0 ) lying above m , the profile of the M 1 -component of the corresponding wave solution is heteroclinic with M 1 (c(−∞)) < M 1 (c(∞)), whereas the one lying below m , the corresponding profile of the M 1 -component is heteroclinic with M 1 (c(−∞)) > M 1 (c(∞)). The aforementioned discussions are illustrated in Figs. 5 and 6.
In contrast to the profile of the M 1 component of wave solutions, the profile of the M 2 component is always monotone increasing, as can be verified by (21) and (49). Further, by adding (18) and (21), we have that Therefore, the profile of the total buffers in calcium-bound forms is always monotone increasing. Finally, let us emphasize that similar monotonicity/nonmonotonicity properties are shared by the profiles of traveling wave solutions to system (5), as it has been shown in Figs. 5 and 6 .

Conclusions
In this work, we considered the properties of traveling wave solutions in systems of reaction-diffusion equations describing the dynamics of calcium ions in the presence of buffer molecules with two binding sites, which can reciprocally influence each other depending on their state, i.e., on whether they are free or occupied. Moreover, in the model we assumed that the binding takes place sequentially-the second after the first. Likewise the unbinding takes place in the inverse direction. It is obvious that, if the binding and unbinding processes at the two sites (in a given buffer molecule) are independent, then the action of 'two-site' buffers coincides with the action of two subpopulations of 'one-site' buffers with halved total densities. However, in general, one site and two site models give different quantitative as well as qualitative results. The proposed model with two site buffering molecules is described by system (5) (Buffered bistable system) or the corresponding system of odes for the traveling wave solutions (10). In Sects. 3.2-3.4, we derive from system (5), using the rapid buffering approximation, the Rapid buffered bistable system (23), (18), (21) and its parabolic counterpart (25), (18), (21). In Sect. 4, we formulate the condition for the existence of traveling waves (with positive speed) independently of how large is the buffer diffusion coefficient (D M ) and the total buffer concentration b 0 . Basing on the asymptotic approximation of the proposed model, we can analyze the action of two-site buffers. In particular, we can compare it with the effect of one-site buffering molecules. In our study we chose the equilibrium level of cytosolic calcium ions (ground state) as c b = 0.2 (μM), but all the conclusions derived in the paper hold qualitatively for all positive values of c b . First of all, as it was mentioned in Sect. 5.2.1, two-site buffers can stop calcium traveling waves for smaller values of the parameter a ∈ (0, 1/2), i.e., for higher excitability, with respect to one-site buffers (0.254 vs. 0.331). In the same subsection we analyze the shape of admissible regions in the K 0 K 1 -space for two chosen excitability parameters a equal to 0.38 (high excitability) and a = 0.43 (low excitability). These shapes can have a very complicated form as it is seen from the middle and right panel of Fig. 3. In particular, there emerges an additional region guaranteeing the advancing wave propagation independently of the value D M b 0 .
In Sect. 5.2.3 we consider the problem of monotonicity of traveling wave profiles as solutions to the asymptotic equations (23), (18) and (21). It thus follows from relation (18) that for certain regions of parameter pairs (K 0 , K 1 ) the profiles of M 1 can be nonmonotone or even monotonically decreasing. Moreover, in Figs. 5, 6 we show numerically that the same monotonicity properties hold also for M 1 profiles of traveling wave solutions to the initial system (5). This consistency can be treated as a numerical proof of the validity of the asymptotic reduction. Finally, two-site buffers can also facilitate the propagation of calcium traveling waves, i.e., induce their propagation, even in the case when it is impossible without their presence. This is shown numerically in Fig. 4 both for the rapid buffered bistable system and for the initial system (5) (Buffered bistable system). This is a highly nonintuitive and unexpected result. Although it is well known that calcium buffers (including the major calcium buffers calmodulin and calretinin) generally have multiple binding sites that are not independent (Starovasnik et al. 1992;Schwaller 2010;Prins and Michalak 2011), the possible effects of this on the speed and existence of physiological waves remain almost entirely unexplored and unknown. Unfortunately, it is difficult to perform an experiment in which the interactions of the individual calcium buffer binding sites can be modulated while leaving the total amount of buffering unchanged; to our knowledge, this has never been attempted. Moreover, such an experiment would have to be performed ensuring that the buffering power remains within a relatively narrow region where the effects on waves of binding site dependence can be seen.
In future work, we intend to generalize our analysis to buffers in which the binding sites not only influence each other, but also have different binding and unbinding coefficients. This however will necessitate considering a system with more equations. and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: Existence of Traveling Waves for Sufficiently Fast Binding and Unbinding Reactions
In this appendix, we will consider the existence of solutions to the system i.e., system (10). In the proof, we will use the modified method applied to a system describing one kind of buffer molecules with one binding site Kaźmierczak and Peradzynski (2011). As it was done in Sect. 3.1, let us rescale the kinetic coefficients of the considered system: Hence, system (A.1) can be written in the form In fact, we will assume that the parameter L is sufficiently large. Let us define From (A.4) and (A.5) we obtain: and, in turn, the expressions for the derivatives of M 1 (·) and M 2 (·). Having the expressions for M 1 (·), M 2 (·), M 1 (·) and M 2 (·), we can write the second and third equation of system (A.1) in the form whereas˜ 1 and˜ 2 are functions of c, c , η 1 , η 2 and v. Finally where θ 1 and θ 2 are defined in (18) and (21), respectively. In (A.8) the functions˜ 1 ,˜ 2 contain the terms v M 1 and v M 2 , hence the terms proportional to c (·), η 1 (·) and η 2 (·) (obtained via the differentiation of (A.6) and (A.7)), as well as the products of these derivatives contained in the functions D M M 1 and D M M 2 .
Next, let us note that c (·) can be obtained by means of the first equation of system (A.1). This equation can be written as

Equation (A.8) can thus be written as
Consequently, we obtain the system: Let us consider the system corresponding to (A.8) for the variables η := η 1 + η 2 , η 2 . It is seen that this system has the form: Let us note that the matrix A d is obtained from the matrix K = A * sec 1 + χφ 1 2χφ 2 χφ 1 1 + 2χφ 2 via the relations We will first analyze the properties of the matrix where φ 10 and φ 20 are defined in (A.10). This matrix can be written as Putting the expressions for φ 10 and φ 20 from (A.10) we obtain: Lemma A.1 There exists a constant c b ∈ (0, 1), such that for all sets of the coefficients κ 0 ± > 0, κ 1 ± > 0, the trace of the matrix B positive.
Proof As it follows from straightforward calculations, the trace of B is equal to: The first component of the sum is obviously positive. To prove the positivity of the second component, let us note that the expression in the main bracket is positive for κ 0 − ≥ 0 and c > 0, unless the coefficient (κ 0 + ) 2 + (−1 + c)κ 0 + κ 1 + + 2κ 1 − κ 1 + is negative. Next, the minimal value of the second-order polynomial with respect to κ 0 Suppose that c < 1. Then, (−1 + c)κ 0 + κ 1 + < 0, but as we have assumed that [(κ 0 The second term of the above sum is positive for c ∈ (0.1716, 5.8284) hence for c ∈ (0.1716, 1). If c ≥ 1, then [(κ 0

Proof By straightforward calculation it follows that
As the eigenvalues of the matrix B are equal to This expression is positive unless the coefficient (−κ 1 − + c(κ 0 + + κ 1 + )) is negative. Then, the minimal value of the quadratic polynomial with respect to κ 0 − is equal to The last expression is positive, because κ 1 − − cκ 0 + > cκ 1 + > 0 for c > 0. It follows that the eigenvalues of A are positive for some d ∈ [0, d 1 ) for some d 1 > 0 sufficiently small. Now, let us consider the case of d · b 0 sufficiently large . Then, by (A.20), the Det(B) grows proportionally to db 0 , whereas (T race(B)) 2 grows proportionally to (db 0 ) 2 . It follows that condition (A.21) is satisfied, so the second part of the thesis of the lemma holds. The lemma is proved. Now, we can generalize the result concerning the eigenvalues of the matrix B expressed in Lemma A.2, to the matrix A d , which in contrast to the matrix B, is dependent on the functions η, η 2 . In this case to retain the positivity of the eigenvalues of A d , we must assume that the C 0 norms of the functions η and η 2 are sufficiently small.
with s 0 sufficiently small (but independent of the parameter L). Then, there exist D 1 > 0 sufficiently small and D 2 > D 1 sufficiently large, such that for db 0 ∈ [0, D 1 ) ∪ (D 2 , ∞), the eigenvalues of A d (c,η,η 2 ) are both positive for all ξ ∈ IR.
Remark For simplicity, in our denotation we have not emphasized in Definition 1 the fact that the eigenvalues depend on the functions c(·),η(·) andη 2 (·). So in fact we should use the notation: Let us start from the analysis of the homogeneous system adjoint to the homogeneous version of system (A.16), i.e., We will show that Eq. (A.22) has no bounded smooth solutions in C 2 (IR) except for zero solutions. Consider the matrix where C 1 and C 2 are the column eigenvectors corresponding to the positive eigenvalues λ 1,2 (ξ ) of the matrix A d (c,η,η 2 )(ξ ), i.e., Remark As the eigenvalues of the matrix A d (c,η,η 2 ), the entries of the matrix S(ξ ), i.e., depend on the functions c,η, andη 2 . So, S(ξ ) = S c,η,η 2 (ξ ).
To be more precise where the matrices Z (ξ ) and Z 1 (ξ ) do not depend on the parameter L. It follows that system (A.28) can be written in the form: (A.29) As system (A.29) is linear, then we can normalize it by the demand that one of the norms: ζ 1 C 2 (IR) or ζ 2 C 2 (IR) is equal to 1 and the other is not bigger than 1. To fix our attention, we will assume that ζ 1 C 1 (IR) = 1 and ζ 2 C 1 (IR) ≤ 1. (A.30) Let us denote Suppose that the function ζ 1 (·) attains a positive maximum at a point ξ * ∈ IR. Then, as Z and Z 1 do not depend on L, ζ 1 (ξ * ) = 0 and ζ 1 (ξ * ) < 0, hence thus, due to the assumed normalization Likewise, if ζ 1 attains a negative minimum at ξ * , then where C is independent of L.
In the same way, we conclude that Remark In analyzing the extrema of the functions ζ 1 and ζ 2 , we assumed that they are attained at finite values of ξ . The same reasoning can be applied also to extrema attained at infinities. Consider, e.g., a positive maximum of ζ 1 achieved at +∞. It can happen, if there exists a sequence of increasing maxima, so it can be analyzed as above. In the other, monotone case, for sufficiently large ξ we must have ζ 1 (ξ ) > 0, ζ 1 (ξ ) → 0 for ξ → ∞, and ζ 1 (ξ ) < 0, so we can also use the arguments applied to the case of finite ξ * . The other cases can be analyzed similarly.
Next, it follows from the assumed normalization that the second derivatives of the functions ζ 1 and ζ 2 are bounded from above and below by numbers independent of L.
Lemma A.6 Consider the homogeneous system corresponding to Eq. (A.22): Suppose that the functions c(·),η andη 2 satisfy the conditions of Lemma A.5. Then, for L sufficiently large, there is no C 2 -bounded nonzero solution to system (A.35).
Now, we will consider the system corresponding to system (A.16): with the functions c(·),η(·) andη 2 fixed and satisfying conditions listed in Lemma A.5. Let us recall that, under these conditions, the matrix A d (c,η,η 2 )(ξ ) can be diagonalized according to (A.24) by means of the matrix S −1 (ξ ) given by (A.26). Thus, multiplying system (A.36) by the matrix S −1 (ξ ), we obtain by means of (A.24) the system: Let us note that the matrices and 1 correspond to matrices Z and Z 1 defined after (A.28) by replacing J (ξ ) by S −1 (ξ ). In this way, by denoting system (A.37) can be written in the form: where θ i j are the entries of the matrix .
To proceed, let us note that system (A.38) can be written by introducing additional dependent variables corresponding to the first derivatives (ν + 1 := ζ and ν 2 := ζ 2 ) as a first order system (of four equations), i.e., System (A.39) can be formally treated as a linear system. For this representation of system (A.16), we can use Lemma 4.2 in Palmer (1984). According to Lemma A.7, the matrix multiplyingL for every x ∈ IR, the 4×4 matrix at the left-hand side of the above system has two positive and two negative eigenvalues (as we assume that the functions c,η andη 2 satisfy the conditions of Lemma A.5), so the dimensions of the stable and unstable subspaces are equal to 2. It follows that the left-hand side operator in (A.39) is Fredholm of index 0. Obviously, Lemma A.5 implies that the system adjoint to the homogeneous counterpart of system (A.39) has no nonzero solutions bounded in C 1 (IR) class. It follows that, for A f and ( 1 , 2 ) T treated as given functions, system (A.39) has a unique C 1 (IR) solution.
Consequently, given A f , 1 and 2 , system (A.38) has a unique C 2 -bounded solution. Let us establish estimates for the derivatives of this solution. To do this we will make stronger assumptions concerning the functionsc,η andη 2 . Namely, we will additionally assume that c C 3 (IR) , η C 3 (IR) and η 2 C 3 (IR) are bounded. (A.40) First, let us estimate first the C 0 -norms of the solutions to system (A.38). Suppose that at ξ = ξ * the function ζ 1 (ξ ) attains a positive maximum. According to the Remark after inequality (A.34) we can confine ourselves to the case ξ * ∈ IR. It follows from the first equation of the system that ζ(ξ * ) ≤ CL −1 . Proceeding in this way, we can prove that ζ 1 (·) C 0 (IR) ≤CL −1 , ζ 2 (·) C 0 (IR) ≤CL −1 .
Moreover, according to Lemma A.5 the above estimates should be written as Next by differentiation of the equations in system (A.38) with respect to ξ , and using the estimates of C 0 -norms, we can conclude that By means of the above estimates, it follows from Eq. (A.38) that ζ 1 (ξ ) and ζ 2 (ξ ) can be bounded from above and below by a constant independent ofL. Consequently, it follows from the differentiation of Eq. (A.38) with respect to ξ that ζ 1 (ξ ) and ζ 2 (ξ ) are bounded uniformly for all ξ ∈ IR. (A.43)

Remark
The differentiation of the equations of system (A.38) is well determined as the entries of the matrix is of C 1 class and the entries of the matrix 1 are of C 2 class, as we assume (see A.40) that the functionsc,η andη 2 are of C 3 class. Next, the eigenvalues λ 1 and λ 2 depend only on these functions (but not of its derivatives) and 1 and 2 are of C 2 class.
In view of the fact that the transformation (η, η 2 ) → (ζ 1 , ζ 2 ) acting from C 2 (IR) to C 2 (IR) is well defined and invertible, then, according to Lemma A.5, the following lemma holds.
Proof The second claim of the lemma can be proved by repeating the reasoning given after system (A.39).
Next, due to the properties of the transformation (η, η 2 ) → (ζ 1 , ζ 2 ), it follows that according to (A.41) and (A.42) we can write Differentiating Eq. (A.16) twice with respect to ξ we obtain Using (A.44), we conclude that the right-hand side of the above equation can be estimated by a constant not depending on L for all ξ ∈ IR. Multiplying the both sides of the last equation by S −1 (ξ ) we obtain the equation: In view of the first remark after (A.42), the second and the third derivatives of the functions η and η 2 are of the order of O(1), hence by denoting We can write Eq. (A.46) as: where the functions T and T 2 are of the order of O(1) and independent of L. It thus follows from the maximum principle and the Remark after (A.34) that Next, using the C 0 (IR) boundedness of the functions η and η 2 and using the arguments from the proof of Lemma A.4, we can show that We have thus proved the following lemma.
Lemma A.9 Suppose that η, η 2 denote a solution to Eq. (A.16), i.e., Suppose that the solution satisfies the estimate η(·) C 2 (IR) ≤ C, η 2 (·) C 2 (IR) ≤ C, where the constant C does not depend on L. Then, the following estimates hold: where the constant C depends on the C 2 properties of the entries of the matrix A d , but does not depend on L as L → ∞.
Let us note that the first equation of system (A.1) can be written as and consequently in the form independent of L: Differentiating (A.6) and (A.7) we can subsequently write it in the form: and φ 1 and φ 2 are defined in (A.10), whereas d M in (A.9). As it follows from (24) and (A.10), where ν 2 is a linear function of η 1 and η 2 , vanishing for η 1 = η 2 = 0. Its form can be deduced from (A.10) to (A.9). Exactly, the same properties have the functions ν 1 and ν 0 . Moreover, as it is easy to check that similar properties are shared by the function c with respect to η 1 , η 2 as well as with respect to their derivatives. Our aim is to prove the existence of a traveling wave solution to system (A.1) for sufficiently large values of the parameter L and to show that as L → ∞ this traveling wave solution tends together with its speed of propagation to the heteroclinic pair given by the asymptotic Eq. (23). Our tool will be the implicit function theorem in Banach spaces of differentiable functions.
Definition 2 For = i = 0, 1, 2, let B i denote the sub-space of functions u(z) belonging to BC i (IR) and tending to finite limits as z → ±∞ together with their derivatives (which tend to zero). Let B i0 denote the subspace of B i consisting of functions u satisfying the condition: The norms in the spaces B j are taken to be In these norms B i and B i0 are Banach spaces.

Proof
The proof of the lemma follows by a straightforward checking.
We have shown the validity of the following statement.
In view of Lemmata A.10-A.12, we can use the implicit function theorem to prove the existence of traveling waves for system (A.1). It thus suffices to show that the linear system D Q| λ=0 [δc, δη, δη 2 , δv] = h c ∈ B 1 δη = h η ∈ B 3 δη 2 = h η 2 ∈ B 3 (A.56) defines an isomorphism between the spaces B 30 × B 3 × B 3 × IR, i.e., (via Theorem 4.2-H in Taylor (1958)) and that the above system is boundedly invertible. Replacing δη by h η and δη 2 by h η 2 in the first equation of the above system becomes Let us recall that according to Lemma 1, w (ξ ) > 0 for all ξ ∈ IR, so due the monotonicity of transformation (27), the profile C(·) of the traveling front to Eq. (A.49) with η ≡ 0, η 2 ≡ 0 (and λ = 0), is monotonically increasing, i.e., C (ξ ) > 0. Moreover, C (ξ ) → 0 as ξ → ±∞. It follows that the homogeneous version of Eq. (A.57) has a unique bounded solution (up to a translation and a multiplicative constant) equal to C (·), because the other linearly independent solution diverges at infinities. However, this function does not belong to the space B 30 , because it does not satisfy condition (A.51). It may be easily shown that a unique (up to a multiplicative constant) solution to a conjugated equation has the form C (ξ ) exp( In view of the fact that C (ξ ) > 0 for all ξ ∈ IR, the last equation is uniquely solvable with respect to δv. Leading order approximations As the functions and 2 are independent of the parameter L, then, according to Lemma A.11 in the first nonzero approximations: