Part IV: Electrophysiology & Signaling | Chapter 3

Calcium Signaling & Dynamics

Calcium buffering, IP$_3$ receptor kinetics, calcium waves, and CICR oscillations — the universal second messenger decoded through biophysics

Why Calcium?

Calcium ions (Ca$^{2+}$) serve as the most versatile intracellular signal in biology. The resting cytosolic concentration ($\sim 100$ nM) is maintained 10,000-fold below the extracellular level ($\sim 1.5$ mM) and 1000-fold below the endoplasmic reticulum (ER) store ($\sim 100$–$500$ $\mu$M). This enormous gradient makes Ca$^{2+}$ an ideal signaling molecule: small fluxes produce large concentration changes.

This chapter derives the quantitative framework for calcium signaling: buffering kinetics that shape transients, IP$_3$ receptor models that control ER release, reaction-diffusion equations for calcium waves, and the bifurcation analysis of calcium oscillations.

1. Calcium Buffering

Most Ca$^{2+}$ entering the cytoplasm is immediately captured by buffer proteins (calmodulin, parvalbumin, calbindin). Only 1–5% of total calcium remains free. Understanding buffering is essential for interpreting every calcium measurement.

Derivation: Rapid Buffer Approximation

Consider a single buffer species B that binds Ca$^{2+}$ with on-rate$k_{\text{on}}$ and off-rate $k_{\text{off}}$:

$$\text{Ca}^{2+} + \text{B} \underset{k_{\text{off}}}{\overset{k_{\text{on}}}{\rightleftharpoons}} \text{CaB}$$

The dissociation constant is $K_d = k_{\text{off}}/k_{\text{on}}$. The total calcium is partitioned:

$$[\text{Ca}]_{\text{total}} = [\text{Ca}^{2+}] + [\text{CaB}]$$

The total buffer is conserved: $[B]_{\text{total}} = [B] + [\text{CaB}]$. At equilibrium:

$$[\text{CaB}] = \frac{[B]_{\text{total}} \cdot [\text{Ca}^{2+}]}{K_d + [\text{Ca}^{2+}]}$$

In the rapid buffer approximation (RBA), we assume the buffer equilibrates much faster than other calcium dynamics. Then for a change$\delta[\text{Ca}]_{\text{total}}$, the fraction appearing as free calcium is:

$$\frac{\delta[\text{Ca}^{2+}]}{\delta[\text{Ca}]_{\text{total}}} = \frac{1}{1 + \kappa}$$

where $\kappa$ is the buffer capacity:

$$\boxed{\kappa = \frac{d[\text{CaB}]}{d[\text{Ca}^{2+}]} = \frac{[B]_{\text{total}} K_d}{(K_d + [\text{Ca}^{2+}])^2}}$$

At resting calcium ($[\text{Ca}^{2+}] \ll K_d$ for most buffers):

$$\kappa \approx \frac{[B]_{\text{total}}}{K_d}$$

Typical values: endogenous buffers give $\kappa \sim 20$–$200$, meaning only 0.5–5% of entering calcium remains free.

For multiple buffers, the total capacity adds:

$$\kappa_{\text{total}} = \sum_j \frac{[B_j]_{\text{total}} K_{d,j}}{(K_{d,j} + [\text{Ca}^{2+}])^2}$$

Derivation: Effective Diffusion Coefficient

Calcium diffusion in the cytoplasm is slowed by buffering. The reaction-diffusion equations are:

$$\frac{\partial [\text{Ca}^{2+}]}{\partial t} = D_{\text{Ca}} \nabla^2 [\text{Ca}^{2+}] - k_{\text{on}} [\text{Ca}^{2+}][B] + k_{\text{off}} [\text{CaB}]$$

$$\frac{\partial [\text{CaB}]}{\partial t} = D_B \nabla^2 [\text{CaB}] + k_{\text{on}} [\text{Ca}^{2+}][B] - k_{\text{off}} [\text{CaB}]$$

In the rapid buffer approximation with an immobile buffer ($D_B = 0$), we combine these equations. Since $[\text{CaB}]$ tracks $[\text{Ca}^{2+}]$ instantaneously:

$$\nabla^2[\text{CaB}] = \kappa \nabla^2[\text{Ca}^{2+}]$$

Adding the two equations:

$$(1 + \kappa) \frac{\partial [\text{Ca}^{2+}]}{\partial t} = (D_{\text{Ca}} + \kappa D_B) \nabla^2 [\text{Ca}^{2+}]$$

The effective diffusion coefficient is:

$$\boxed{D_{\text{eff}} = \frac{D_{\text{Ca}} + \kappa D_B}{1 + \kappa}}$$

For an immobile buffer ($D_B = 0$):

$$D_{\text{eff}} = \frac{D_{\text{Ca}}}{1 + \kappa}$$

With $D_{\text{Ca}} \approx 220$ $\mu$m$^2$/s and$\kappa \approx 100$, the effective diffusion is only$D_{\text{eff}} \approx 2$ $\mu$m$^2$/s — over 100-fold slower than free diffusion. This means calcium signals are highly localised: the characteristic spread distance in time $t$ is$\ell = \sqrt{2 D_{\text{eff}} t}$. For $t = 10$ ms,$\ell \approx 0.2$ $\mu$m. Mobile buffers (like calbindin with$D_B \approx 20$ $\mu$m$^2$/s) increase the effective diffusion and expand the signaling range.

2. IP$_3$ Receptor Kinetics

The IP$_3$ receptor (IP$_3$R) is the primary calcium release channel on the endoplasmic reticulum. Its unique property is biphasic regulation by calcium: low Ca$^{2+}$ activates, high Ca$^{2+}$ inhibits. This bell-shaped calcium dependence is the engine driving calcium oscillations.

Derivation: De Young-Keizer Model

The De Young-Keizer (also called Li-Rinzel) model treats each IP$_3$R subunit as having three independent binding sites:

• IP$_3$ binding site: Fraction bound = $m_{\infty} = [\text{IP}_3]/([\text{IP}_3] + d_1)$
• Ca$^{2+}$ activation site: Fraction bound = $n_{\infty} = [\text{Ca}^{2+}]/([\text{Ca}^{2+}] + d_5)$
• Ca$^{2+}$ inhibition site: Fraction bound = $h$, which evolves dynamically

The IP$_3$R is a tetramer; each subunit must have IP$_3$ bound and Ca$^{2+}$ activating but not Ca$^{2+}$ inhibiting. The open probability for a single subunit is proportional to $m_{\infty} \cdot n_{\infty} \cdot h$, and for the tetramer:

$$\boxed{P_O = \left(\frac{[\text{IP}_3]}{[\text{IP}_3] + d_1} \cdot \frac{[\text{Ca}^{2+}]}{[\text{Ca}^{2+}] + d_5} \cdot h\right)^3}$$

The inhibition variable $h$ evolves slowly:

$$\frac{dh}{dt} = a_2 \left(d_2 \frac{[\text{IP}_3] + d_1}{[\text{IP}_3] + d_3}\right)(1-h) - a_2 [\text{Ca}^{2+}] h$$

At steady state:

$$h_{\infty} = \frac{d_2([\text{IP}_3] + d_1)}{d_2([\text{IP}_3] + d_1) + [\text{Ca}^{2+}]([\text{IP}_3] + d_3)}$$

The time constant of $h$ is:

$$\tau_h = \frac{1}{a_2\left(d_2\frac{[\text{IP}_3]+d_1}{[\text{IP}_3]+d_3} + [\text{Ca}^{2+}]\right)}$$

Bell-shaped Ca$^{2+}$ dependence: The steady-state open probability $P_O^{\infty} \propto n_{\infty}^3 h_{\infty}^3$ first increases with $[\text{Ca}^{2+}]$ (activation via $n_{\infty}$) then decreases (inhibition via $h_{\infty}$). This produces the characteristic bell shape with peak near $[\text{Ca}^{2+}] \approx 0.2$–$0.3$ $\mu$M. Typical parameter values: $d_1 = 0.13$ $\mu$M,$d_2 = 1.049$ $\mu$M, $d_3 = 0.9434$ $\mu$M,$d_5 = 0.08234$ $\mu$M, $a_2 = 0.2$ $\mu$M$^{-1}$s$^{-1}$.

3. Calcium Waves

Calcium waves propagate through cells and tissues via a mechanism of calcium release triggering further release at adjacent sites (calcium-induced calcium release, CICR). These waves transmit information over distances of 10–1000 $\mu$m at speeds of 10–100 $\mu$m/s.

Derivation: Fire-Diffuse-Fire Model

The fire-diffuse-fire (FDF) model treats IP$_3$R clusters as discrete point sources separated by distance $d$. Each cluster "fires" (releases calcium) when the local Ca$^{2+}$ exceeds a threshold $C_{\text{th}}$, then becomes refractory for time $\tau_{\text{ref}}$.

After a cluster fires and releases flux $\sigma$ (amount of Ca$^{2+}$), the calcium diffuses to neighbouring clusters. In 1D, the calcium profile from a point source at $x = 0$ is:

$$C(x, t) = \frac{\sigma}{\sqrt{4\pi D_{\text{eff}} t}} \exp\left(-\frac{x^2}{4 D_{\text{eff}} t}\right)$$

The next cluster at distance $d$ reaches threshold when$C(d, t^*) = C_{\text{th}}$. Solving for $t^*$:

$$C_{\text{th}} = \frac{\sigma}{\sqrt{4\pi D_{\text{eff}} t^*}} \exp\left(-\frac{d^2}{4 D_{\text{eff}} t^*}\right)$$

The wave speed is $v = d/t^*$. For large $\sigma$ (strong release), the dominant balance gives:

$$t^* \approx \frac{d^2}{4 D_{\text{eff}} \ln(\sigma/(C_{\text{th}}\sqrt{4\pi D_{\text{eff}} t^*}))}$$

In the continuum limit ($d \to 0$ with release per unit length held constant), the wave becomes a travelling front in a reaction-diffusion equation. For a simple CICR model with release rate $k_{\text{release}}$:

$$\boxed{v = \sqrt{2 D_{\text{eff}} \cdot k_{\text{release}}}}$$

This Luther-speed formula (analogous to Fisher-KPP wave speed) gives a lower bound on wave velocity.

Condition for propagation failure: The wave fails if the calcium released by one cluster is insufficient to trigger the next. This occurs when:

$$\frac{\sigma}{\sqrt{4\pi D_{\text{eff}} t_{\text{opt}}}} < C_{\text{th}}, \quad t_{\text{opt}} = \frac{d^2}{2 D_{\text{eff}}}$$

$$\boxed{\sigma < C_{\text{th}} d \sqrt{2\pi e} \quad (\text{propagation failure condition})}$$

where $t_{\text{opt}}$ is the time at which the concentration at distance$d$ is maximised (optimising the Gaussian). This shows that propagation failure occurs when clusters are too far apart, release is too weak, or the threshold is too high.

4. Calcium Oscillations

Many cell types exhibit repetitive calcium transients (oscillations) with periods ranging from seconds to minutes. These oscillations arise from the interplay of positive feedback (CICR) and negative feedback (Ca$^{2+}$-dependent inactivation of IP$_3$R, SERCA pump refilling of ER stores).

Derivation: Two-Pool Model (ER + Cytosol)

The minimal model for calcium oscillations tracks two variables: cytosolic Ca$^{2+}$ concentration $c$ and ER Ca$^{2+}$concentration $c_{\text{ER}}$ (or equivalently the inhibition variable $h$). Using the Li-Rinzel reduction of the De Young-Keizer model:

$$\frac{dc}{dt} = J_{\text{rel}} - J_{\text{pump}} + J_{\text{leak}}$$

where the fluxes are:

IP$_3$R release flux:

$$J_{\text{rel}} = k_f \left(\frac{[\text{IP}_3]}{[\text{IP}_3] + d_1}\right)^3 \left(\frac{c}{c + d_5}\right)^3 h^3 (c_{\text{ER}} - c)$$

SERCA pump flux (Hill kinetics):

$$J_{\text{pump}} = V_{\text{SERCA}} \frac{c^2}{K_{\text{SERCA}}^2 + c^2}$$

Passive leak:

$$J_{\text{leak}} = k_{\text{leak}}(c_{\text{ER}} - c)$$

The ER concentration evolves (with volume ratio $\gamma = V_{\text{cyt}}/V_{\text{ER}}$):

$$\frac{dc_{\text{ER}}}{dt} = \gamma(J_{\text{pump}} - J_{\text{rel}} - J_{\text{leak}})$$

For total calcium conservation ($c_{\text{total}} = c + c_{\text{ER}}/\gamma$ = const), we can eliminate $c_{\text{ER}}$ and work with $(c, h)$ as the two dynamical variables. The $h$ equation provides the slow negative feedback:

$$\frac{dh}{dt} = \frac{h_{\infty}(c) - h}{\tau_h(c)}$$

This forms a classic slow-fast system: $c$ is the fast variable (with CICR positive feedback) and $h$ is the slow variable (Ca$^{2+}$-dependent inhibition). The mechanism of oscillation is: (1) low Ca$^{2+}$, $h$is large, IP$_3$R is primed; (2) small trigger releases Ca$^{2+}$ via CICR (positive feedback through $n_{\infty}$); (3) high Ca$^{2+}$drives $h$ down (inhibition); (4) IP$_3$R closes, SERCA pumps Ca$^{2+}$ back to ER; (5) low Ca$^{2+}$ allows $h$ to recover, cycle repeats.

Hopf Bifurcation Analysis

The two-variable $(c, h)$ system undergoes a Hopf bifurcation as the IP$_3$concentration increases. At the fixed point $(c^*, h^*)$, the Jacobian is:

$$J = \begin{pmatrix} \partial \dot{c}/\partial c & \partial \dot{c}/\partial h \\ \partial \dot{h}/\partial c & \partial \dot{h}/\partial h \end{pmatrix}$$

The key entries are:

• $\partial \dot{c}/\partial c$: Contains the CICR positive feedback (positive contribution from $\partial J_{\text{rel}}/\partial c > 0$) minus the SERCA pump removal (negative)
• $\partial \dot{c}/\partial h > 0$: More release when inhibition gate is open
• $\partial \dot{h}/\partial c < 0$: Higher calcium causes more inhibition (the negative feedback)
• $\partial \dot{h}/\partial h < 0$: The inhibition variable relaxes to its steady state

The Hopf bifurcation occurs when the trace of $J$ crosses zero:

$$\text{tr}(J) = \frac{\partial \dot{c}}{\partial c} + \frac{\partial \dot{h}}{\partial h} = 0$$

Since $\partial \dot{h}/\partial h < 0$ (always stabilising), the bifurcation occurs when $\partial \dot{c}/\partial c$ becomes sufficiently positive, i.e., when the CICR positive feedback becomes strong enough to overcome the pump and leak stabilisation.

Increasing [IP$_3$] strengthens the release flux (more IP$_3$R are activated), making $\partial \dot{c}/\partial c$ more positive. At a critical [IP$_3$]$^*$, the system transitions from a stable fixed point to a limit cycle (oscillations). At even higher [IP$_3$], a second Hopf bifurcation returns the system to a stable (but elevated) steady state, as the ER is depleted and CICR positive feedback weakens.

5. Applications

Muscle Contraction

In skeletal muscle, the ryanodine receptor (RyR) mediates CICR from the sarcoplasmic reticulum. The Ca$^{2+}$ transient activates troponin C, enabling actin-myosin cross-bridge cycling. The Hill model of force-calcium relationship:

$$F = F_{\max} \frac{[\text{Ca}^{2+}]^{n_H}}{K_{0.5}^{n_H} + [\text{Ca}^{2+}]^{n_H}}$$

with $n_H \approx 2.5$ and $K_{0.5} \approx 0.6$ $\mu$M, provides a steep, switch-like activation. The steep Hill coefficient arises from cooperative calcium binding to the troponin complex.

Fertilisation

Sperm-egg fusion triggers a dramatic calcium wave that sweeps across the egg, preventing polyspermy and initiating embryonic development. In mammalian eggs, the sperm introduces PLC$\zeta$, which generates IP$_3$, triggering repetitive calcium oscillations lasting several hours. The frequency encoding of calcium signals activates CaMKII, which drives cell cycle progression.

Cardiac Arrhythmias

Abnormal calcium handling underlies many cardiac arrhythmias. Spontaneous calcium release from the SR (calcium sparks that propagate as waves) activates the Na$^+$/Ca$^{2+}$exchanger (NCX), generating a depolarising transient inward current ($I_{\text{ti}}$) that can trigger delayed after-depolarisations (DADs). When DADs reach threshold, they produce triggered beats that initiate arrhythmias.

Neuronal Plasticity

The direction and magnitude of synaptic plasticity depends on postsynaptic calcium dynamics. The BCM (Bienenstock-Cooper-Munro) theory proposes a modification threshold$\theta_m$ that itself depends on calcium history:

• Modest Ca$^{2+}$ elevation ($\sim 0.3$–$0.5$ $\mu$M) activates calcineurin, triggering long-term depression (LTD)
• Large Ca$^{2+}$ elevation ($> 1$ $\mu$M) preferentially activates CaMKII, triggering long-term potentiation (LTP)
• The different calcium thresholds arise from different $K_d$ values: calcineurin has higher affinity ($K_d \sim 0.1$ $\mu$M) than CaMKII ($K_d \sim 1$ $\mu$M)

6. Python Simulations

Calcium Oscillations & Buffering Analysis

Python

script.py209 lines

import numpy as np

# =========================================================
# CALCIUM OSCILLATION SIMULATION
# Li-Rinzel reduction of De Young-Keizer model
# =========================================================

print("=" * 60)
print("CALCIUM OSCILLATIONS: Li-Rinzel Model")
print("=" * 60)

# Parameters (Li-Rinzel model)
# IP3R parameters
d1 = 0.13       # IP3 dissociation constant (uM)
d2 = 1.049      # Ca inhibition dissociation (uM)
d3 = 0.9434     # IP3 dissociation (inhibition site, uM)
d5 = 0.08234    # Ca activation dissociation (uM)
a2 = 0.2        # Ca inhibition rate (uM^-1 s^-1)

# Flux parameters
kf = 1.11       # max IP3R flux rate (s^-1)
V_SERCA = 0.9   # max SERCA pump rate (uM/s)
K_SERCA = 0.1   # SERCA half-activation (uM)
k_leak = 0.002  # passive leak rate (s^-1)

# Cell geometry
gamma = 5.405   # cytoplasm-to-ER volume ratio
c_total = 2.0   # total Ca (uM, conserved in cytosol + ER)

print("\nModel parameters:")
print(f"  d1={d1}, d2={d2}, d3={d3}, d5={d5} uM")
print(f"  a2={a2} uM^-1 s^-1")
print(f"  kf={kf} s^-1, V_SERCA={V_SERCA} uM/s, K_SERCA={K_SERCA} uM")
print(f"  k_leak={k_leak} s^-1, gamma={gamma}")

def m_inf(IP3):
    return IP3 / (IP3 + d1)

def n_inf(c):
    return c / (c + d5)

def h_inf(c, IP3):
    num = d2 * (IP3 + d1)
    den = d2 * (IP3 + d1) + c * (IP3 + d3)
    return num / den

def tau_h(c, IP3):
    return 1.0 / (a2 * (d2 * (IP3 + d1) / (IP3 + d3) + c))

def J_release(c, h, IP3, c_ER):
    m = m_inf(IP3)
    n = n_inf(c)
    Po = (m * n * h)**3
    return kf * Po * (c_ER - c)

def J_pump(c):
    return V_SERCA * c**2 / (K_SERCA**2 + c**2)

def J_leak_fn(c, c_ER):
    return k_leak * (c_ER - c)

# --- Simulate oscillations at several IP3 levels ---
print("\n" + "=" * 60)
print("OSCILLATION DYNAMICS AT DIFFERENT [IP3]")
print("=" * 60)

dt = 0.005  # seconds
T = 120.0   # seconds
N = int(T / dt)
t = np.arange(N) * dt

IP3_values = [0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 3.0]

for IP3 in IP3_values:
    c = np.zeros(N)
    h = np.zeros(N)

# Initial conditions
    c[0] = 0.05  # resting cytosolic Ca (uM)
    h[0] = h_inf(c[0], IP3)

for k in range(N - 1):
        c_ER = (c_total - c[k]) * gamma  # from conservation

J_rel = J_release(c[k], h[k], IP3, c_ER)
        J_p = J_pump(c[k])
        J_l = J_leak_fn(c[k], c_ER)

dc = (J_rel - J_p + J_l) * dt
        dh = (h_inf(c[k], IP3) - h[k]) / tau_h(c[k], IP3) * dt

c[k+1] = max(c[k] + dc, 0.001)
        h[k+1] = min(max(h[k] + dh, 0.0), 1.0)

# Analyse: find peaks in second half
    c_second = c[N//2:]
    t_second = t[N//2:]

# Detect peaks (simple method)
    peaks = []
    for k in range(1, len(c_second) - 1):
        if c_second[k] > c_second[k-1] and c_second[k] > c_second[k+1] and c_second[k] > 0.15:
            peaks.append(k)

# Remove peaks too close together
    if len(peaks) > 1:
        filtered = [peaks[0]]
        for p in peaks[1:]:
            if (p - filtered[-1]) * dt > 1.0:
                filtered.append(p)
        peaks = filtered

if len(peaks) >= 2:
        periods = np.diff(np.array(peaks) * dt)
        period = np.mean(periods)
        amplitude = np.max(c_second) - np.min(c_second)
        print(f"  [IP3]={IP3:.1f} uM: OSCILLATING, period={period:.1f} s, "
              f"amp={amplitude:.3f} uM, peak={np.max(c_second):.3f} uM")
    else:
        print(f"  [IP3]={IP3:.1f} uM: STEADY STATE, c_ss={np.mean(c_second[-1000:]):.4f} uM")

# --- Detailed simulation at IP3 = 0.6 uM ---
print("\n" + "=" * 60)
print("DETAILED OSCILLATION at [IP3] = 0.6 uM")
print("=" * 60)

IP3_det = 0.6
c_det = np.zeros(N)
h_det = np.zeros(N)
c_ER_trace = np.zeros(N)

c_det[0] = 0.05
h_det[0] = h_inf(c_det[0], IP3_det)

for k in range(N - 1):
    c_ER_val = (c_total - c_det[k]) * gamma
    c_ER_trace[k] = c_ER_val

J_rel = J_release(c_det[k], h_det[k], IP3_det, c_ER_val)
    J_p = J_pump(c_det[k])
    J_l = J_leak_fn(c_det[k], c_ER_val)

dc = (J_rel - J_p + J_l) * dt
    dh = (h_inf(c_det[k], IP3_det) - h_det[k]) / tau_h(c_det[k], IP3_det) * dt

c_det[k+1] = max(c_det[k] + dc, 0.001)
    h_det[k+1] = min(max(h_det[k] + dh, 0.0), 1.0)

c_ER_trace[-1] = (c_total - c_det[-1]) * gamma

print(f"\nTime course (selected points):")
print(f"  {'t (s)':>8} {'[Ca] (uM)':>10} {'h':>8} {'[Ca]_ER':>10}")
for idx in range(0, N, N // 20):
    print(f"  {t[idx]:>8.1f} {c_det[idx]:>10.4f} {h_det[idx]:>8.4f} {c_ER_trace[idx]:>10.2f}")

# IP3R open probability vs Ca at steady state
print("\n" + "=" * 60)
print("IP3R OPEN PROBABILITY: Bell-shaped Ca dependence")
print("=" * 60)

ca_range = np.logspace(-2, 1, 50)
print(f"\n{'[Ca] (uM)':>12} {'P_open (IP3=0.3)':>18} {'P_open (IP3=1.0)':>18}")
print("-" * 52)

for ca_val in [0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 1.0, 2.0, 5.0, 10.0]:
    for ip3_val in [0.3, 1.0]:
        m = m_inf(ip3_val)
        n = n_inf(ca_val)
        h_ss = h_inf(ca_val, ip3_val)
        Po = (m * n * h_ss)**3
        if ip3_val == 0.3:
            print(f"  {ca_val:>10.3f}", end="")
            print(f"  {Po:>16.6f}", end="")
        else:
            print(f"  {Po:>16.6f}")

# --- Buffering analysis ---
print("\n" + "=" * 60)
print("CALCIUM BUFFERING ANALYSIS")
print("=" * 60)

# Buffer parameters
buffers = [
    ("Calmodulin", 24.0, 1.0, 0),       # B_total (uM), Kd (uM), D_B (um^2/s)
    ("Parvalbumin", 50.0, 0.05, 0),      # slow, high affinity
    ("Calbindin-D28k", 150.0, 0.5, 20),  # mobile
    ("BAPTA (exo)", 100.0, 0.22, 95),    # fast, mobile chelator
    ("EGTA (exo)", 100.0, 0.15, 113),    # slow chelator
]

D_Ca = 220.0  # um^2/s

print(f"\n{'Buffer':>16} {'[B] (uM)':>10} {'Kd (uM)':>8} {'kappa':>8} {'D_eff (um^2/s)':>15}")
print("-" * 62)

ca_rest = 0.1  # uM

for name, B_tot, Kd, D_B in buffers:
    kappa = B_tot * Kd / (Kd + ca_rest)**2
    D_eff = (D_Ca + kappa * D_B) / (1 + kappa)
    print(f"  {name:>14} {B_tot:>10.1f} {Kd:>8.2f} {kappa:>8.1f} {D_eff:>15.2f}")

total_kappa = sum(B * K / (K + ca_rest)**2 for _, B, K, _ in buffers[:3])
D_eff_total = D_Ca / (1 + total_kappa)  # assuming immobile
print(f"\n  Endogenous buffer capacity (CaM+PV+CB): kappa = {total_kappa:.1f}")
print(f"  D_eff (immobile buffers only) = {D_eff_total:.2f} um^2/s")
print(f"  Fraction of Ca free = 1/(1+kappa) = {1/(1+total_kappa)*100:.2f}%")
print(f"  Spread distance in 10 ms = sqrt(2*D_eff*t) = {np.sqrt(2*D_eff_total*0.01)*1000:.0f} nm")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Calcium Wave Propagation & Bifurcation Analysis

Python

script.py273 lines

import numpy as np

# =========================================================
# CALCIUM WAVE PROPAGATION: Fire-Diffuse-Fire Model
# =========================================================

print("=" * 60)
print("CALCIUM WAVE: Fire-Diffuse-Fire Model")
print("=" * 60)

# Parameters
D_eff = 5.0       # effective diffusion coefficient (um^2/s)
d_cluster = 2.0   # distance between IP3R clusters (um)
C_th = 0.3        # threshold for triggering release (uM)
sigma = 1.0       # amount of Ca released per cluster (uM * um in 1D)
tau_release = 0.1 # duration of release event (s)
C_rest = 0.05     # resting Ca concentration (uM)

print(f"\nParameters:")
print(f"  D_eff = {D_eff} um^2/s")
print(f"  Cluster spacing d = {d_cluster} um")
print(f"  Threshold C_th = {C_th} uM")
print(f"  Release amount sigma = {sigma} uM*um")

# Analytical wave speed estimates
# Luther formula for continuum limit
k_release_eff = sigma / (d_cluster * tau_release)  # effective release rate
v_Luther = np.sqrt(2 * D_eff * k_release_eff)
print(f"\nContinuum (Luther) wave speed: v = sqrt(2*D*k) = {v_Luther:.1f} um/s")

# FDF optimal arrival time
t_opt = d_cluster**2 / (2 * D_eff)
C_max = sigma / np.sqrt(4 * np.pi * D_eff * t_opt)
print(f"\nFDF analysis:")
print(f"  Optimal arrival time t* = d^2/(2D) = {t_opt:.4f} s")
print(f"  Max concentration at neighbor = {C_max:.4f} uM")
print(f"  Threshold = {C_th} uM")
if C_max > C_th:
    print(f"  Wave PROPAGATES (C_max > C_th)")
else:
    print(f"  Wave FAILS (C_max < C_th)")

# Propagation failure condition
sigma_crit = C_th * d_cluster * np.sqrt(2 * np.pi * np.e)
print(f"  Critical sigma for propagation = {sigma_crit:.4f} uM*um")
print(f"  Actual sigma = {sigma:.4f} uM*um")

# --- Simulate 1D wave propagation ---
print("\n" + "=" * 60)
print("1D WAVE SIMULATION")
print("=" * 60)

# Spatial grid
L = 50.0      # total length (um)
dx = 0.1      # spatial step (um)
Nx = int(L / dx)
x = np.arange(Nx) * dx

# Time parameters
dt_sim = 0.0002  # time step (s)
T_sim = 3.0      # total time (s)
Nt = int(T_sim / dt_sim)

# Initialize
c = np.ones(Nx) * C_rest
c_prev = c.copy()

# Simple CICR reaction-diffusion model
# dc/dt = D * d2c/dx2 + J_release(c) - J_pump(c)

# Release: Hill function with threshold-like behavior (CICR)
def release_rate(c_local, c_ER_local):
    # IP3R-like release: positive feedback
    n_hill = 4
    K_act = 0.3   # uM, activation threshold
    V_rel = 5.0   # max release rate (uM/s)
    activation = c_local**n_hill / (K_act**n_hill + c_local**n_hill)
    driving = max(c_ER_local - c_local, 0)
    return V_rel * activation * driving

def pump_rate(c_local):
    V_pump = 1.0  # uM/s
    K_pump = 0.2  # uM
    return V_pump * c_local**2 / (K_pump**2 + c_local**2)

# ER concentration (simplified: uniform initially)
c_ER = 5.0  # uM (ER store)
c_ER_arr = np.ones(Nx) * c_ER

# Track cluster firing
cluster_positions = np.arange(0, L, d_cluster)
cluster_indices = (cluster_positions / dx).astype(int)
cluster_fired = np.zeros(len(cluster_indices), dtype=bool)
cluster_fire_time = np.full(len(cluster_indices), -1.0)

# Trigger first cluster
c[0:int(1.0/dx)] = 1.0  # initial stimulus

# Store snapshots
snapshot_times = [0, 0.2, 0.5, 1.0, 1.5, 2.0, 2.5]
snapshots = {}

alpha_diff = D_eff * dt_sim / dx**2

print(f"\nGrid: {Nx} points, dx={dx} um")
print(f"Time: {Nt} steps, dt={dt_sim*1000:.2f} ms")
print(f"Diffusion stability number: {alpha_diff:.4f} (should be < 0.5)")

for step in range(Nt):
    t_now = step * dt_sim

# Store snapshots
    for st in snapshot_times:
        if abs(t_now - st) < dt_sim / 2:
            snapshots[st] = c.copy()

# Diffusion (explicit finite differences)
    c_new = c.copy()
    for i in range(1, Nx - 1):
        diffusion = alpha_diff * (c[i+1] - 2*c[i] + c[i-1])
        # Reaction
        J_rel = release_rate(c[i], c_ER_arr[i])
        J_p = pump_rate(c[i])
        reaction = (J_rel - J_p) * dt_sim

c_new[i] = c[i] + diffusion + reaction
        c_new[i] = max(c_new[i], 0.001)

# Update ER (simplified)
        c_ER_arr[i] = max(c_ER_arr[i] - J_rel * dt_sim * 0.2, 0.1)

# Boundary conditions (no flux)
    c_new[0] = c_new[1]
    c_new[-1] = c_new[-2]
    c = c_new

# Print snapshots
print("\nCalcium wave snapshots:")
for st in sorted(snapshots.keys()):
    snap = snapshots[st]
    # Find wave front (where c > 0.5 * max)
    wave_front = -1
    c_max = np.max(snap)
    for i in range(Nx - 1, 0, -1):
        if snap[i] > 0.3:
            wave_front = x[i]
            break

print(f"\n  t = {st:.1f} s: peak [Ca] = {c_max:.3f} uM, wave front at x = {wave_front:.1f} um")
    # Print profile at selected x positions
    print(f"    {'x (um)':>8} {'[Ca] (uM)':>10}")
    for xi in range(0, Nx, Nx // 10):
        if snap[xi] > C_rest * 1.1:
            print(f"    {x[xi]:>8.1f} {snap[xi]:>10.4f}")

# Estimate wave speed from front positions
front_positions = []
front_times = []
for st in sorted(snapshots.keys()):
    if st > 0:
        snap = snapshots[st]
        for i in range(Nx - 1, 0, -1):
            if snap[i] > 0.3:
                front_positions.append(x[i])
                front_times.append(st)
                break

if len(front_positions) >= 2:
    speeds = []
    for i in range(1, len(front_positions)):
        dx_front = front_positions[i] - front_positions[i-1]
        dt_front = front_times[i] - front_times[i-1]
        if dt_front > 0:
            speeds.append(dx_front / dt_front)

if speeds:
        mean_speed = np.mean(speeds)
        print(f"\nEstimated wave speed: {mean_speed:.1f} um/s")
        print(f"Luther formula prediction: {v_Luther:.1f} um/s")

# --- Bifurcation: IP3 vs oscillation amplitude ---
print("\n" + "=" * 60)
print("BIFURCATION: IP3 CONCENTRATION vs OSCILLATION")
print("=" * 60)

# Use simplified Li-Rinzel model
d1 = 0.13
d2 = 1.049
d3 = 0.9434
d5 = 0.08234
a2_lr = 0.2
kf_lr = 1.11
V_S = 0.9
K_S = 0.1
k_l = 0.002
gam = 5.405
c_tot = 2.0

IP3_scan = np.linspace(0.1, 3.0, 50)
amp_results = []
period_results = []

dt_bif = 0.01
T_bif = 200.0
N_bif = int(T_bif / dt_bif)

for ip3 in IP3_scan:
    c_b = 0.05
    h_b = d2 * (ip3 + d1) / (d2 * (ip3 + d1) + c_b * (ip3 + d3))

c_trace = np.zeros(N_bif)
    for k in range(N_bif):
        c_trace[k] = c_b
        c_ER_b = (c_tot - c_b) * gam

m_b = ip3 / (ip3 + d1)
        n_b = c_b / (c_b + d5)
        Po_b = (m_b * n_b * h_b)**3
        J_r = kf_lr * Po_b * (c_ER_b - c_b)
        J_p = V_S * c_b**2 / (K_S**2 + c_b**2)
        J_l = k_l * (c_ER_b - c_b)

h_ss = d2 * (ip3 + d1) / (d2 * (ip3 + d1) + c_b * (ip3 + d3))
        th = 1.0 / (a2_lr * (d2 * (ip3 + d1) / (ip3 + d3) + c_b))

c_b = max(c_b + (J_r - J_p + J_l) * dt_bif, 0.001)
        h_b = h_b + (h_ss - h_b) / th * dt_bif
        h_b = min(max(h_b, 0), 1)

# Analyse last half
    c_last = c_trace[N_bif//2:]
    c_min = np.min(c_last)
    c_max_val = np.max(c_last)
    amp = c_max_val - c_min
    amp_results.append(amp)

# Detect period
    peaks_b = []
    for k in range(1, len(c_last) - 1):
        if c_last[k] > c_last[k-1] and c_last[k] > c_last[k+1] and c_last[k] > c_min + 0.1:
            peaks_b.append(k)
    if len(peaks_b) >= 2:
        filtered_b = [peaks_b[0]]
        for p in peaks_b[1:]:
            if (p - filtered_b[-1]) * dt_bif > 0.5:
                filtered_b.append(p)
        if len(filtered_b) >= 2:
            per = np.mean(np.diff(filtered_b)) * dt_bif
            period_results.append(per)
        else:
            period_results.append(0)
    else:
        period_results.append(0)

print(f"\n{'[IP3] (uM)':>12} {'Amplitude':>10} {'Period (s)':>12} {'Regime':>12}")
print("-" * 50)
for k in range(0, len(IP3_scan), 4):
    ip3 = IP3_scan[k]
    amp = amp_results[k]
    per = period_results[k]
    regime = "oscillating" if amp > 0.05 else "stable"
    per_str = f"{per:.1f}" if per > 0 else "---"
    print(f"  {ip3:>10.2f} {amp:>10.4f} {per_str:>12} {regime:>12}")

print("\n" + "=" * 60)
print("SUMMARY")
print("=" * 60)
print("1. Calcium buffering: kappa ~ 20-200, D_eff = D_Ca/(1+kappa)")
print("2. IP3R has bell-shaped Ca dependence (activation + inhibition)")
print("3. Calcium waves: v ~ sqrt(2*D*k), fail if spacing too large")
print("4. Oscillations arise from CICR + slow Ca-dependent inhibition")
print("5. Hopf bifurcation: oscillations appear in window of [IP3]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Chapter Summary

• Calcium buffering: Buffer capacity $\kappa = [B]K_d/(K_d+[\text{Ca}])^2$ determines that only $1/(1+\kappa)$ of calcium remains free. Effective diffusion $D_{\text{eff}} = D_{\text{Ca}}/(1+\kappa)$ for immobile buffers.
• IP$_3$ receptor: The De Young-Keizer model captures the bell-shaped Ca$^{2+}$ dependence through three binding sites (IP$_3$, activating Ca$^{2+}$, inhibiting Ca$^{2+}$). The slow inhibition variable $h$ drives oscillations.
• Calcium waves: The fire-diffuse-fire model gives wave speed $v \sim \sqrt{2D_{\text{eff}} k_{\text{release}}}$. Propagation fails when cluster spacing exceeds $\sigma/(C_{\text{th}}\sqrt{2\pi e})$.
• Oscillations: CICR positive feedback and slow Ca$^{2+}$-dependent inhibition create a Hopf bifurcation: oscillations exist in a window of IP$_3$ concentrations.
• Applications: Muscle contraction, fertilisation waves, cardiac arrhythmias, and synaptic plasticity all depend on precise calcium dynamics.

← Neural Biophysics Mechanotransduction →