Long-range diffusion of bound Ca 2+: time-dependent WKBJ solutions in one- and two-dimensional imaging

1 Introduction

The Michaelis and Menten (1913) equation and its associated quasi-steady-state approximation (Segel and Slemrod, 1989; Bersani et al., 2015; Choi et al., 2017) have served for more than a century as fundamental tools for describing biochemical and biophysical reactions. Despite numerous refinements, these formulations generally describe either spatially homogeneous systems or highly localized reactions in which diffusion is neglected.

For intracellular Ca²⁺ signaling, Wagner and Keizer (1994) incorporated Michaelis–Menten kinetics, also known as the rapid buffer approximation (RBA) for Ca²⁺ buffering, into a reaction–diffusion (RD) framework and derived a non-linear equation for total Ca²⁺ concentration, which they analyzed numerically and under several limiting approximations. Mironova and Mironov (2008) found that total Ca²⁺ can be well described using the logarithmic approximation and derived analytical formulas for transient free Ca²⁺ dynamics. They also emphasized that Ca²⁺ is transported predominantly in its bound form within the cytoplasm.

Many endogenous Ca²⁺ buffers diffuse slowly or are effectively immobilized (Zhou and Neher, 1993; Schwaller, 2010). Mironov (2023a) analyzed the case of irreversible Ca²⁺ binding and derived analytical expressions that predicted a localized accumulation of free Ca²⁺ around a single channel. Ca²⁺-bound buffer species were found to propagate away from it much faster. Since endogenous buffers and sensors are usually large proteins with limited mobility, the propagation of bound Ca²⁺ may represent an efficient mechanism for long-range intracellular signaling. In this view, Ca²⁺-bound molecules, rather than free Ca²⁺ itself, may constitute the principal carriers of information between distant intracellular sites.

At first sight, the hypothesis may appear inconsistent with the widely accepted concept of stationary Ca²⁺ nanodomains surrounding single channels (Neher, 1986; Eisner et al., 2023). However, these perspectives are not mutually exclusive, as discussed below.

In the present study, we revisit Michaelis–Menten-type kinetics by explicitly incorporating diffusion and derive analytical solutions using the Liouville–Green (LG) transformation, commonly known in quantum mechanics as the Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) approximation. Although WKBJ methods are routinely applied to stationary problems, their use in non-stationary RD systems has remained limited (Krause et al., 2020).

To maintain the accessibility of the main text, all mathematical derivations are presented in Sections 5.1 and 5.2. We demonstrate that the leading-order term of the LG/WKBJ expansion accurately reproduces multidimensional diffusion with and without reactions and provides simple analytical expressions for Ca²⁺ buffering. Experimental measurements obtained in one-dimensional single-channel recordings and two-dimensional neuronal preparations support the theoretical predictions.

Overall, our findings suggest that Ca²⁺-bound species can propagate more rapidly and over greater distances than localized free Ca²⁺ signals, potentially influencing the timing and spatial organization of intracellular signaling. In addition to Ca²⁺ dynamics, the analytical framework developed here offers a straightforward approach for incorporating diffusion into Michaelis–Menten-type reaction systems.

2 Methods

2.1 Ethical approval

All animal experiments were approved by the Lower Saxony State Office for Consumer Protection and Food Safety and conducted in accordance with institutional and national guidelines.

2.2 Cell preparation

Cultured hippocampal neurons and chick dorsal root ganglion (DRG) neurons were prepared as previously described (Mironov, 1990; Mironov et al., 1993, respectively). DRG neurons were used 6–12 h after plating, when they had fully adhered to the substrate and exhibited a flattened circular morphology with a thickness of approximately 2–3 μm. At this time in culture, spontaneous synaptic activity was absent, minimizing contributions from neuronal network activity.

2.3 Solutions

In single-channel measurements, the solutions contained 30 mM Tris buffer (pH 7.4) together with either 154 mM NaCl (pipette) or 88 mM CaCl₂ (bath). Fluo-4 and Fluo-4 dextran (30 kDa) were used at final concentrations of 0.2 mM. In whole-cell recordings, the intrapipette solution contained CsMeSO4 (92 mM), CsCl (43 mM), TEA-Cl (5 mM), MgCl2 (1 mM), HEPES (10 mM), ATP (4 mM), and GTP (0.4 mM); a pH of 7.4 (titrated with CsOH) was used. The indicators were added at 200 μM.

2.4 Patch-clamp recording and imaging

Electrophysiological recordings and fluorescence imaging were performed as previously described. α-Synuclein channels were incorporated into inside-out membrane patches excised from the cultured hippocampal neurons (Mironov, 2015). Opening of individual channels generated outward Ca²⁺ currents that produced localized increases in Ca²⁺-bound fluorescence within the recording pipette.

In whole-cell recordings, voltage steps were applied from −60 to −10 mV and evoked Ca²⁺ currents that were predominantly carried by L-type Ca²⁺ channels.

Fluorescence imaging was performed using total internal reflection fluorescence microscopy (TIRF) combined with high-speed electron-multiplying charge-coupled device (EMCCD) acquisition. In whole-cell recordings, the experimental protocol—control, BayK 8044, and nifedipine (see Methods)—was applied when the L-type current showed no significant rundown within the first 5 min.

For imaging, an upright microscope (Axioscope 2, Zeiss) with a 63 × objective lens (numerical aperture, 1.4) was used. A laser beam delivered 488-nm light from below through a prism (TIRF Labs, Inc., Cary, NC) at an angle appropriate for TIRF. The 525-nm emission passed through a dichroic mirror was filtered at 535 nm. Images were captured using a cooled iXON Ultra 888 EMCCD camera (ANDOR, Oxford Instruments) operated with ANDOR software. In the one-dimensional (1D) experiments, line-scan mode (512 × 2 pixels at 12-bit resolution) was used, corresponding to an acquisition time of 0.1 ms in Crop mode. In whole-cell recordings, the recorded area was 128 × 128 pixels, and the acquisition time was 1.5 ms. Fluo-4 fluorescence increases >30-fold after Ca²⁺ binding that minimized out-of-focus effects. The lateral resolution, measured using 0.1-μm TetraSpeck beads (ThermoFisher, Germar, Germany), yielded a lateral point-spread function of 0.18 μm. The TIRF depth was approximately 200 nm, allowing the flattened DRG cells in culture to be considered planar objects and the results obtained for Ca²⁺ diffusion in circular geometry to be applied. Image analysis was performed offline using MetaMorph software. Images were background-subtracted and median-filtered before analysis. Fluorescence changes are reported as relative changes (ΔF/F₀).

2.5 Statistics

Each experiment was repeated in at least four cells obtained from different preparations. Data are presented as mean ± SD. Statistical significance was assessed using paired t-tests where appropriate, with significance accepted at p < 0.05.

2.6 Reaction–diffusion model

Ca²⁺ diffusion in the presence of a single buffer was described by coupled reaction–diffusion Equation 1:Governing free Ca²⁺, free buffer (F), and Ca²⁺-bound buffer (B) concentrations. Here, the Laplacian is written in n-dimensional form (Carslaw and Jaeger, 1986; Crank, 1979). The model incorporates diffusion coefficients and reversible binding kinetics characterized by association (kon) and dissociation (koff) rate constants representative of cytoplasmic Ca²⁺-binding proteins. The Ca²⁺-buffering constants (Table 1) have typical values for cytoplasmic calcium-binding proteins: kon = 3^.108 M⁻¹ s⁻¹ corresponds to diffusion-limited binding, Bo = 0.2 mM, and DCa = 220 μm²/s. Numerical values for the model parameters are summarized in Table 1.

[ Ca ] t = D Ca ( [ Ca ] rr + n − 1 r [ Ca ] r ) − k on [ Ca ] [ F ] + k off [ CaB ] (1)

[ F ] t = D buffer ( [ F ] rr + n − 1 r [ F ] r ) − k on [ Ca ] [ F ] + k off [ CaB ]

Table 1. Basic definitions and model parameters.

Analytical solutions were derived using the LG/WKBJ formalism (Sections 5.1 and 5.2). To leading order, the resulting expressions provide compact closed-form approximations for free Ca²⁺ and bound buffer distributions in both one-dimensional and multidimensional geometries. For the 1D case:Where A is the concentration of Ca²⁺ at the boundary; z = x/2√t is a self-similarity (Boltzmann) variable, Q = √(z2 + 1), and D = √(1 - d)/2, with the normalized parameter d = Dbuffer/DCa. The concentration of free buffer is defined byAnd:Where the non-dimensional dissociation constant K = Kd/Bd= 0.003 (Table 1). The expressions are consistent with the RBA and Michaelis–Menten formalism. For radial diffusion (the 3D case), all concentrations are simply divided by r; for the 2D case, a numerical solution for free Ca²⁺ is preferable because the analytical formulas are cumbersome. The data in Section 5.3 show that Equations 2–4 reproduce well results obtained by direct numerical solution of RD-PDE reaction–diffusion equations.

[ C a 2 + ] = A exp ( − zQ − z 2 / 2 ) √ ( z + Q ) (2)

[ F ] = B o K ( [ Ca 2 + ] + K ) (3)

[ CaB ] = B o [ Ca ] ( [ Ca ] + K ) (4)

3 Results

3.1 Theoretical prediction of one- and two-dimensional Ca²⁺ profiles

Fluorescence measurements report Ca²⁺ bound to the indicator rather than free Ca²⁺ itself, and the theoretical analysis therefore focuses on the dynamics of the bound species. The LG/WKBJ formalism was used to derive analytical expressions for the spatiotemporal distributions of free Ca²⁺ and Ca²⁺-bound buffer species in Section 5.2. Figure 1 compares solutions obtained for slowly diffusing buffers (relative diffusion coefficient d = Dbuffer//DCa = 0.03) and more rapidly diffusing buffers (d = 0.8; e.g., DFluo−4 = 180 μm2/s, Mironova and Mironov, 2008). In both cases, the leading-order WKBJ approximation closely matches the solution including the first-order correction, demonstrating rapid convergence of the expansion and indicating that the leading term alone can provide a good approximation.

For slowly diffusing buffers, free Ca²⁺ spreads progressively from the source as the buffer becomes locally depleted. In contrast, when buffer mobility is closer to that of free Ca²⁺, both free and bound Ca²⁺ distributions remain substantially more confined. This behavior is consistent with the narrow Ca²⁺ nanodomains typically inferred from measurements obtained with small organic fluorescent indicators, whose diffusion coefficients are only moderately lower than that of free Ca²⁺. Comparison of one-dimensional and circular geometries shows that the predicted free Ca²⁺ profiles resemble each other at distances greater than 1 μm and times exceeding 2 ms. Deviations near the source arise from the geometry-dependent Laplacian term and become negligible at larger distances.

For a sustained Ca²⁺ source, free Ca²⁺ concentrations increase continuously because progressive buffer saturation reduces local buffering capacity. In contrast, the distributions of Ca²⁺-bound and free buffer evolve as traveling fronts that propagate away from the source. Calculated velocities are similar for linear and circular geometries (Table 2), indicating that this behavior is largely independent of spatial dimensionality.

Table 2. Estimated velocities (μm/s) of Ca²⁺-bound buffer waves.

Theoretical velocities from front position at 12 ms from data shown in Figures 2, 3. In the second column the numbers are obtained for basic (default) model parameters. In next columns the velocities are obtained in calculation runs, when of only one parameter in the basic set was modified as indicated in the first line. Experimental velocity estimates give the mean wave speeds ± SD as measured for Ca2+-bound buffer (organic indicator dye, fluo-4(X)) in 1D- (single α-Syn channel, Figure 5) and 2D-systems (DRG neurons, Figure 6). Experimental estimates were obtained from FWHM (midpoints) of linescans in Figures 5, 6 divided by elapsed time. The data were analyzed for four different patches with α-Syn channels and four whole-cell patched cells. Control values vs. BAYK data had P values equal to 0.0006 (Fluo-4) and 0.0088 (Fluo-4-dextran) respectively.

The calculations shown in Figure 1 assume sustained Ca²⁺ flux. For channels with shorter open times, such as voltage-gated Ca²⁺ channels, the expansion of profiles would be terminated upon channel closure and therefore acquire spatially narrower patterns. Under these conditions, the apparent width of the free Ca²⁺ domain is determined jointly by diffusion and channel open time, whereas the distribution of bound Ca²⁺ additionally reflects the kinetics of buffer binding. Consequently, the duration of Ca²⁺ influx may critically influence downstream signaling by determining whether sufficient sensor activation can occur before channel closure.

One possible test of the proposed model would involve examining the correlation between the opening duration of Ca²⁺ channels and the activity of Ca²⁺-dependent channels when they are present in the same patch. A previous study (Mironov, 2015) found that synuclein channel activity did not alter the activity of delayed rectifier channels (Kdr) but modulated ATP-sensitive K⁺ (KATP) channel activity through changes in submembrane ion homeostasis. It would be interesting to correlate the opening of Ca²⁺ channels and Ca²⁺-activated K⁺ (KCa) channels. However, a preparation in which KCa channels are abundant would be required; this is not the case for the hippocampal neurons used in this study. The long-term effects of bound Ca²+ could be expected for non-selective Ca²+-permeable channels as NMDA and transient receptor potential (TRP) channels (Mironov, 2008).

Figure 1 illustrates the effects of key model parameters. After normalization to the Ca²⁺ concentration at the source, the shape of the free Ca²⁺ profile is largely independent of influx amplitude. In contrast, increasing the source Ca²⁺ concentration markedly accelerates propagation of the bound species while preserving similar wave morphology in both one- and two-dimensional geometries. This behavior arises because elevated Ca²⁺ rapidly saturates local buffer molecules, causing the depletion front to propagate outward.

The dissociation constant likewise influences wave propagation. Lower K_d values produce faster-moving fronts, whereas larger K_d values yield slower and more spatially confined propagation owing to more rapid equilibration between free and bound states. However, across all parameter combinations examined, the qualitative propagating behavior of Ca²⁺-bound species remains robust.

The sigmoidal shape of the propagating fronts follows directly from the analytical solution. Under the rapid-buffer approximation, the concentration of bound buffer can be expressed as a Michaelis–Menten-type function of free Ca²⁺, Equation 4. The resulting bound fraction naturally assumes a hyperbolic tangent–like profile, with its midpoint determined primarily by the effective dissociation constant. Under the RBA (Michaelis and Menten, 1913) approximation, the concentration of Ca²⁺-bound buffer can be expressed as:

[ CaB ] B o = e ( − s − lnK ) ( e ( − s − lnK ) + 1 ) = [ tanh ( ( s + lnK ) / 2 ) + 1 ] (5)

Thus, the logarithm of the dissociation constant serves as the midpoint of the wave of bound Ca²⁺. When K increases, its logarithm becomes more negative, which leads to denser spacing of the profiles (e.g. Figure 2C vs. Figure 2D).

3.2 One-dimensional propagation of Ca²⁺-bound species generated by single α-synuclein channels

To test the theoretical predictions experimentally, one-dimensional Ca²⁺ profiles were generated using outward currents through individual α-synuclein (α-Syn) channels incorporated into excised inside-out membrane patches. The channels are approximately equally permeable to K⁺, Na⁺, and Ca²⁺ (Mironov, 2015), but because the bath contains excess Ca²⁺, the recorded ‘outward’ current represents an almost pure Ca²⁺ current through a single α-Syn channel. This preparation offers the advantage that channel openings frequently persist for several milliseconds, permitting direct observation of the temporal evolution of Ca²⁺-dependent fluorescence signals.

The patches were held at a potential that generated a pure 0.2-pA Ca²⁺ current, matching a typical single L-type channel current (Mironov and Richter, 1998) and corresponding to [Ca2⁺]_o/Bo = 0.1 for the default parameter set (Table 1). Figure 4 shows that opening of a single channel produced localized increases in Fluo-4 or Fluo-4-dextran fluorescence corresponding to the formation of Ca²⁺-bound indicator molecules. Kymographs revealed a progressive expansion of the fluorescence signal from the pipette tip toward the pipette interior. Line scans extracted at successive time points demonstrated that the fluorescence front propagated continuously away from the source rather than remaining spatially stationary.

The observed propagation is consistent with the predicted wave behavior of Ca²⁺-bound species. Although fluorescence decayed rapidly, following channel closure because of diffusion and ongoing binding reactions, the spatial evolution during channel opening was well described by a propagating front. Mean experimental velocities are summarized in Table 2, and these velocities agree qualitatively with theoretical predictions.

3.3 Two-dimensional profiles of bound Ca²⁺ in the neuronal soma

The theoretical predictions were further evaluated in a two-dimensional cellular geometry using cultured dorsal root ganglion neurons loaded intracellularly with Fluo-4 or Fluo-4-dextran. Whole-cell Ca²⁺ currents were evoked by depolarizing voltage steps, producing reproducible increases in intracellular fluorescence.

At the beginning of the experiment, the voltage pulses produced an increase in the middle of the cell due to excitation of the basal membrane. When the cell was ‘trained’ with 100-ms-long pulses delivered at 0.1 Hz to deplete interstitial space, this offset was markedly reduced, and the lateral transients originating at the cell borders were seen more clearly. Immediately after cell excitation, fluorescence first increased beneath the plasma membrane before progressively spreading into the bulk cytoplasm. Line-scan analysis demonstrated that the resulting Ca²⁺-bound fluorescence propagated inward as a traveling front (Figures 3,4).

Pharmacological manipulation of Ca²⁺ influx produced the changes predicted by the analytical model. Enhancement of L-type Ca²⁺ channel activity with Bay K 8644 increased current amplitude and accelerated propagation of the fluorescence front, whereas inhibition with nifedipine reduced both current amplitude and wave velocity. These effects were observed for both indicator dyes but were particularly pronounced with the slowly diffusing dextran-conjugated probe. Front position changed in proportion to √t, which is characteristic of a diffusion process (the self-similarity variable is z = x/2√t; see Equation 2 above and Section 5.1.).

The experimentally determined propagation velocities (Table 2) followed the theoretical trends predicted by the WKBJ analysis: larger Ca²⁺ influxes produced faster-moving waves, whereas reduced influx slowed propagation. Moreover, the overall similarity between the one-dimensional and two-dimensional measurements supports the conclusion that the spread of Ca²⁺-bound species represents a robust consequence of RD dynamics rather than a feature specific to a particular experimental geometry.

4 Discussion

The present study complements and extends two widely used concepts in biochemistry and biophysics: the Michaelis and Menten (1913) formalism and the Excess Buffer Approximation (EBA) introduced by Neher (1986). The Michaelis–Menten approach describes local substrate–enzyme interactions but does not account for the diffusion of either the substrate or the enzyme. Similarly, the EBA assumes that Ca²⁺ buffers are present in excess and therefore remain far from saturation. This assumption permits linearization of the original RD system (Equation 1), reducing it to second-order ordinary differential equation. The resulting solution predicts the formation of localized Ca²⁺ nanodomains around persistently open Ca²⁺ channels. This concept is frequently invoked to explain a variety of physiological phenomena, as recently reviewed by Eisner et al. (2023).

The majority of theoretical studies of intracellular Ca²⁺ signaling have focused on exponentially decaying steady-state Ca²⁺ distributions (Neher, 1986; Naraghi and Neher, 1997; Smith et al., 2001). However, these analyses critically depend on the assumption of excess buffer. When this constraint is removed, the analytical solution for irreversible Ca²⁺ binding with nearly equal Ca²⁺ and buffer diffusion coefficients predicts the appearance of periodic Ca²⁺ patterns under specific conditions (Mironov, 2019). Exact solutions were also found for an immobile buffer that binds Ca²⁺ irreversibly (Mironov, 2023a; the RBA does not apply in this case). An important conclusion emerging from these analyses is that free Ca²⁺ and Ca²⁺-bound species may exhibit markedly different dynamic behaviors. While free Ca²⁺ accumulates relatively slowly around an individual channel, Ca²⁺ bound to a buffer or sensor can propagate away from the source at substantially higher velocities (Mironov, 2023a). Consequently, the effective signal carriers—the Ca²⁺-bound buffer or sensor molecules—may travel both faster and farther than the signaling ion itself.

To fully understand the information conveyed by intracellular Ca²⁺ signals, it is necessary to consider not only the distribution of free Ca²⁺ but also the generation and propagation of Ca²⁺-bound buffers and sensors. Fluorescence measurements report Ca²⁺ bound to the indicator rather than free Ca²⁺ itself, and the theoretical analysis therefore focused on the dynamics of the bound species. In the theoretical part of this study, we used the RBA, which is identical to Michaelis–Menten kinetics widely used in biochemistry. Wagner and Keizer (1994) were the first to apply the RBA to intracellular Ca²⁺ buffering and derived a non-linear RD equation for total cellular Ca²⁺. Their analysis relies primarily on numerical solutions and limiting cases. However, approximate analytical solutions are possible under the RBA when a logarithmic approximation of total Ca²⁺ is used (Mironova and Mironov, 2008).

The present study further extends this concept by explicitly incorporating diffusion into Michaelis–Menten-type reaction kinetics. To derive analytical solutions, the LG transformation was used, as described in Sections 5.1 and 5.2. In quantum mechanics, this approach is widely known as the WKBJ method and is commonly used to solve the stationary Schrödinger equation. Applications of WKBJ techniques to time-dependent RD systems are comparatively rare. Among the few examples are the studies by Krause et al. (2020), although these authors used a stationary WKBJ formalism in which temporal dependence entered only indirectly.

To avoid overburdening the main text, all mathematical derivations and numerical analyses are presented in Sections 5.1 and 5.2. As an initial validation of the approach, the WKBJ formalism was applied to multidimensional diffusion and to diffusion coupled with a first-order reaction corresponding to the time-dependent version of the EBA (Section 5.1). The resulting solutions accurately reproduced established analytical results. Remarkably, the leading-order WKBJ approximation alone proved sufficient to achieve excellent agreement with exact solutions (Section 5.1; Figure 5).

The same methodology was subsequently extended to the Ca²⁺–buffer system (Section 5.2). The resulting derivations yield simple and transparent analytical expressions, as summarized in Equations 2–4. Once again, the leading-order WKBJ term appeared sufficient to reproduce the behavior of multidimensional RD systems (Figure 1). Although only a single buffer species was considered here, the formalism can readily be extended to multiple buffers with different mobilities, as previously proposed for related RBA models (Mironova and Mironov, 2008).

The theoretical analysis identifies buffer mobility as a major determinant of intracellular Ca²⁺ patterns (Figure 1). Contrary to the predictions of the classical EBA, these patterns are not stationary, even under conditions of constant Ca²⁺ influx. Owing to progressive buffer saturation, free Ca²⁺ distributions continuously expand away from the channel source in both linear (Figures 1A,B) and circular (Figures 1C,D) geometries. Even more striking behavior is observed for Ca²⁺-bound species, whose distributions spread over substantially larger distances.

For buffer diffusion coefficients approaching that of free Ca²⁺ (Figure 1B), the distributions of both free and bound Ca²⁺ remain spatially confined, thereby supporting the classical concept of Ca²⁺ nanodomains (Neher, 1986). This paradigm was used to interpret measurements performed using synthetic fluorescent indicator dyes whose diffusion coefficients are close to that of free Ca²⁺. Endogenous cytoplasmic Ca²⁺ buffers such as calmodulin are larger molecules and possess substantially lower diffusion coefficients; consequently, they are often treated as effectively immobile buffers (Zhou and Neher, 1993; Schwaller, 2010). Both the theoretical predictions presented here (Figures 1, 2) and the corresponding experimental observations suggest that, under conditions of low buffer mobility, Ca²⁺-bound buffer and sensor molecules are not confined to localized domains but instead propagate progressively throughout the cytoplasm.

5 A novel mathematical treatment of RD systems using the WKBJ formalism

5.1 Time-dependent solutions of simple diffusion equations

Before applying the LG/WKBJ formalism to Ca²⁺ buffering, we first validate the approach using diffusion and RD equations for which exact analytical solutions are available. This comparison demonstrates that the leading-order WKBJ approximation accurately reproduces the known analytical solutions and provides the foundation for the subsequent analysis.

Consider a non-stationary (time-dependent) case of simple diffusion:Where subscripts t and r denote the time and spatial derivatives, respectively, and Δ represents the n-dimensional Laplacian:

c t = Δc (A1)

Δc = c rr + n − 1 r c r

Used in Equation 1 (Crank, 1979; Carslaw and Jaeger, 1986). A straightforward way to obtaining the solution is to introduce the Boltzmann (Crank, 1979) or self-similarity variable z = r/2√t (Polyanin and Zaitsev, 2004). Equation A1 then transforms into the following ordinary differential equation (ODE):

c ″ + ( 2 z + n − 1 z c ′ ) = 0 (A2)

Where primes indicate a derivative with respect to z. For further analysis, the equation is transformed into a canonical (LG) form through a change of variable (Polyanin and Zaitsev, 2004):That gives:Where:

c = u exp ( − ∫ 0 z [ z + n − 1 z ] dz )

u ″ − Q ( z ) u = 0 (A3)

Q ( z ) = ( 2 z + ( n − 1 ) z 2 ) 2 / 4 + ( 2 − ( n − 1 ) z 2 ) / 2

Equation A3 is simply a stationary Schrödinger equation with an effective potential Q(z). Following Bender and Orszag (1999), we rewrite it as.

h 2 u ” − Q ( z ) u = 0 (A4)

Where h is a small parameter. The dependent variable is then presented as an exponential, whose argument is expanded in powers of h as

u = exp ( s + hS + … h ) (A5)

Next, substitute (Equation A5) into (Equation A4), divide out the exponential factors, and collect the terms multiplied by the same power of h. The first two equations have solutions (Bender and Orszag, 1999):Where s and S are the first and second terms in the WKBJ expansions, respectively (Equation A8). In most WKBJ applications, further correction terms are small, and only two terms in the expansion of the exponent are sufficient for most purposes (Bender and Orszag, 1999). This also applies in the present study, as shown below. For example, when Q(z) = k, the solution given by Equation A6a,

s = ± ∫ 0 z √ Q dz (A6a)

S = − ln ( Q ) / 4 (A6b)

S = exp ( − kx )

Is identical to the analytical solution of ODE c” = kc, which appears as a steady-state solution under the EBA (Neher, 1986).

For a physically relevant solution, only the exponential with a negative argument (minus sign) in (Equation A6a) is used. For multidimensional diffusion, the equations, potentials Q(z), and leading solutions lnc = s(z) are given in Table 3. In the 1D case, the integral (Equation A6a) is

∫ 0 z √ Q dz = [ zQ / 2 + ln ( z Q + 1 ) / 2 ] (A7)

Table 3. Explicit expressions of functions used in derivation of leading term in WKBJ solutions for multi-dimensional diffusion.

*erfc ( z ) = ∫ z ∞ exp ( − z 2 ) dz, a complementary error function.

For circular diffusion, direct integration of the polynomial Q(z) = √(z2 + 2–1/4z2) is possible, but the result is cumbersome and involves elliptic integrals. Therefore, the original Equation A6 was instead numerically integrated by the method of lines (Berezin and Zhidkov, 1965; Mironov, 2023a), as briefly described in Section 5.3. In the 3D case, the solution is simply the 1D result divided by r.

Figure 4 shows that the leading WKBJ term (Equation A6a) closely reproduces the exact solutions (left panels), and the first-order contribution (Equation A6b) does not modify the results substantially.

Next, consider the simplest RD equation:

c t = △ c − kc (A8)

Compared with pure diffusion, it contains a spatially uniform first-order reaction, whose rate constant is k. Equation A8 corresponds to the time-dependent excess buffer approximation (EBA), which describes irreversible Ca²⁺ binding by the buffer with rate constant k = k_onB_o (Table 1). After transformation into a self-similar form, the equation corresponds to (Equation A3), with 4ktc on the right-hand side. The solutions are readily obtained from (Equation A10) in Table 3 with this term added to Q(z), as determined for pure diffusion. The exact solution (Mironov, 2019) is:

c = [ e − x erfc ( z − t ) + e + x erfc ( z + t ) ] / 2

This solution is again closely reproduced by the leading term (Figure 5, right panels B, D, and F).

5.2 Time-dependent analytical solutions of Ca²⁺ buffering using WKBJ formalism

In all derivations below, we use the normalized equation derived from Equation 1 in Methods. Here, the relative buffer diffusion coefficient is d = Dbuffer/DCa, and the intrinsic time and distance units are τ = 1/konBo = 50 μs and ro = √DCa/Bokon = 105 μm, respectively (Table 1; see also Mironov, 2019, 2023a). The concentrations below are divided by the total buffer concentration Bo = 0.2 mM (Table 1). The system of RD equations for Ca²⁺ buffering isWhere c, f, and b define the concentrations of free Ca²⁺, free buffer, and bound buffer, respectively (all normalized to Bo). Note that b = 1 – f. The normalized dissociation constant is Kd = koff/konBo = 0.003 (Table 1).

c t = Δc − cf + K ( 1 − f ) (B1)

f t = dΔf − c f + K ( 1 − f )

The first two terms describe diffusion, and the last two terms represent the forward and backward steps in the reversible reaction. It is convenient to reorganize the reaction term as:

R = − cf + K ( 1 − f ) = − ( c + K ) f + K = − Cf + K

Where a new Ca²⁺-dependent variable C = c + K is introduced.

In the ‘multiple-scales’ representation (Bender and Orszag, 1999; see also Equation A4 above), we need to solve the system.

h C t = h 2 △ C − Cf + K / h (B2)

h f t = h 2 dΔf − Cf + K / h

When the reaction term:Is expanded in exponential form at leading order h⁻¹, the reaction term R = 0, owing to the identity.

R = − exp ( s h + S + φ h + Φ + … ) + exp ( lnK / h ) (B3)

φ + s = lnK (B4)

Or

f = K / ( K + c ) (B5)

Which describes the equilibrium between Ca²⁺ and the buffer, as presumed in the RBA or Michaelis–Menten formalism.

Expanding the concentrations as exponentials, C = exp(s/h + S + .) and f = exp(φ/h + Φ + .), the system (Equation B3) can be written as:

h s t = h 2 △ s + h s r 2 + R 1 (B6)

h φ t = h 2 d △ φ + h φ r 2 + R 2

Where the reaction terms are:

R 1 = − exp ( φ / h ) + exp ( − s / h + lnK / h )

R 2 = − exp ( s / h ) + exp ( − φ / h + lnK / h )

It follows from (Equation B4) that R1 = R2 = 0. Subtracting the two equations in (Equation B6), we obtain a single equation:

( s − φ ) t = s r 2 − d φ r 2 + hΔ ( s − dφ ) (B7)

It also follows from (Equation B5) that, for both time and spatial derivatives, φ’ = −s’, which simplifies (Equation B8) to:

2 s t = ( 1 − d ) s r 2 (B8)

This is an equation for pure diffusion with an apparent diffusion coefficient D = (1 - d)/2. Its solution for different dimensions can be obtained from line Table 3, with the z variable scaled by √[(1 - d)/2]. Note that the case d = 1 is critical (see also the end of Section 5.3 below) and should be treated separately using, e.g., the full Equation B1, whose solution was obtained previously (Mironov, 2019).

Subtracting the two equations in (Equation B4) at the next order h¹ gives:

( S − Φ ) t = Δ ( s − dφ ) + 2 ( s r S r + D φ r Φ r ) (B9)

To express Φ through S, we write:From which Φ derivatives can be defined as

R = Φ + S = s r 2 − s t

Φ ′ = R ′ − S ′ (B10)

Inserting (Equation B10) into (Equation B9), we obtain a first-order PDE that may be solved analytically (Polyanin et al., 2002), but the explicit solution is rather cumbersome. Therefore, (Equation B9) was solved numerically. The data plotted in Figure 1 in the main text indicate no significant modification of the solution given by the leading WKBJ term (Equation B8).

Explicit formulas for the leading WKBJ term in the 1D formalism are listed in Section 2.6 in Methods, Equations 2–4. The 3D results can be readily obtained from the 1D case by dividing the solution by r. The next section and Figure 6 show that the WKBJ ansatz closely reproduces numerical solutions of the full RD system (Equation B1), indicating that the WKBJ approach offers a simple way to incorporate the effects of diffusion into the Michaelis–Menten (or RBA) formalism. The complex is then transformed irreversibly to the product of reaction showing the progression of the reaction.

5.3 Application of the method of lines to solving RD equations

The algorithm is illustrated for a typical RD equation.

c t = D Δc − f ( c ) (C1)

For diffusion in the presence of a sink, whose concentration dependence may be non-linear. For pure diffusion, the algorithm was described by Berezin and Zhidkov (1965). Below, the algorithm is presented following Mironov (2023b), with appropriate modifications for a non-linear sink. Briefly, write (Equation C1) in discrete form as:

D c n + 2 t − 2 c n + 1 t + c n t d x 2 = [ c n + 1 t − c n + 1 t − Δt dt ] − f ( c n + 1 t − Δt ) (C2)

The solution proceeds in two steps. First, we assume a linear relationship between ctn and ctn₊₁:

c n t = C n t c n + 1 t + ϕ n t (C3)

Setting γ = D dt/dx², we present (Equation C2) as:

γ(c ^t_{n + 2} - 2c ^t_{n + 1} + c^t_n) = (c^t_{n + 1} - c^t-Δt_{n + 1}) - f ( c n + 1 t − Δt ) dt

As in (Equation C3), we can derive a linear relationship between ctn+2 and ctn+1 as:

c t n + 1 = C t n c t n + 2 + φ t n + 1 (C4)

Where:

C t n + 1 = 1 / Den

Den = [ 2 − C t n + 1 / γ ]

φ t n + 1 = φ t n + ( c t − Δt n + 1 − f t − Δ t n + 1 dt ) / Den (C5)

In numerical calculations, the procedure starts from the initial conditions, e.g.:

c ‘ o = 0 , and c o = A

And subsequently evaluates all coefficients C_n and φ^t_n using (Equation C5). The limiting value of cend is then set to 0 or another constant and we move backward to obtain all concentrations c_n at a given t from Equation C3. The two progressions, Equations C4, C5, explicitly define both the forward and backward steps, which explains the name ‘chasing’ for the algorithm (Berezin and Zhidkov, 1965). The algorithm has been tested for many ordinary differential equations (ODEs) and accurately recovered exact solutions (Mironov, 2023b). Compared with other approaches designed to solve RD equations, the main advantage of the method is the ability to use larger time steps without losing accuracy or stability.

The method of lines was used to solve the full Equation B1. Figure 6 presents the results of calculations for the sample cases shown in Figures 1, 2 above. The solutions disagree slightly at d > 0.9, close to the critical value d = 1. In this case, it is better to solve the full system (Mironov, 2019). Nonetheless, the close agreement between the results obtained analytically and numerically indicates that the LG-WKBJ ansatz may be used as a simple analytical alternative to the numerical solution of RD equations.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Luigi Catacuzzeno, University of Perugia, Italy

Reviewed by: Nicolangelo Iannella, University of Oslo, Norway

Tatsat Banerjee, Johns Hopkins University, United States