r"""
Phase boundaries and capture island analysis.
This module implements models for analyzing phase boundaries and capture islands
in the parameter space of the gating system.
Theoretical Foundation
----------------------
The bell-shaped gating function :math:`\mathcal{G}(\phi) = \theta(1-\theta)`
produces **dome-shaped phase boundaries** in concentration space.
For Langmuir isotherm:
.. math::
\mathcal{G}(\phi) = \frac{\phi/K_d}{(1+\phi/K_d)^2}
Given a threshold :math:`g \in (0, 1/4]`, solving :math:`\mathcal{G}(\phi) = g`
yields two solutions:
.. math::
\phi_{\pm}(g) = K_d \cdot \frac{1 - 2g \pm \sqrt{1 - 4g}}{2g}
**Properties:**
- Real solutions exist only for :math:`g \leq 1/4` (peak value)
- Symmetry: :math:`\phi_-(g) \cdot \phi_+(g) = K_d^2` (for :math:`n=1`)
- Peak at :math:`\phi = K_d` where :math:`\mathcal{G} = 1/4`
Capture Island Definition
-------------------------
The "gated regime" is defined by timescale bracketing:
.. math::
\tau_{\text{net}} \lesssim \tau_{\text{esc}}(\phi) \lesssim \tau_{\text{obs}}
This yields concentration intervals forming a closed "band" around
:math:`\phi \approx K_d`, bounded on both low and high concentration sides.
This **re-entrant** behavior is a direct consequence of the bell-shaped
gating function.
References
----------
- Eq. (S63): Gating function definition
- Eq. (S64)-(S66): Boundary solutions for Langmuir
- Eq. (S67): Concentration interval for capture island
"""
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
[docs]@dataclass(frozen=True)
class DomePhaseBoundaries:
r"""
Analytic concentration boundaries for the bell-shaped gating function.
For Langmuir (and Hill via substitution t=(phi/K)^n), the symmetric gating is:
.. math::
\mathcal G(\phi)=\theta(1-\theta)=\frac{t}{(1+t)^2},\quad t=(\phi/K)^n
Given a threshold :math:`g\in(0,1/4]`, solving :math:`\mathcal G=g` yields:
.. math::
t_{\pm}(g)=\frac{1-2g\pm\sqrt{1-4g}}{2g}
hence:
.. math::
\phi_{\pm}(g)=K\,t_{\pm}(g)^{1/n}. \quad (S66\ \text{and Hill extension})
Parameters
----------
K
Dissociation constant K_d (same unit as phi).
n
Hill coefficient. n=1 reduces to Langmuir.
Notes
-----
- Solutions exist only for g <= 1/4 (peak value).
- For n=1, symmetry implies :math:`\phi_-(g)\phi_+(g)=K^2`.
"""
K: float
n: float = 1.0
[docs] def t_solutions(self, g: float) -> tuple[float, float]:
"""
Solve for t in g = t/(1+t)^2.
Parameters
----------
g
Threshold in (0, 1/4].
Returns
-------
t_minus, t_plus
Two positive solutions with t_minus <= 1 <= t_plus.
"""
g = float(g)
if not (0.0 < g <= 0.25):
raise ValueError("g must be in (0, 1/4].")
disc = 1.0 - 4.0 * g
s = np.sqrt(disc)
t_minus = (1.0 - 2.0 * g - s) / (2.0 * g)
t_plus = (1.0 - 2.0 * g + s) / (2.0 * g)
return float(t_minus), float(t_plus)
[docs] def phi_solutions(self, g: float) -> tuple[float, float]:
"""
Solve for phi in G(phi)=g.
Parameters
----------
g
Threshold in (0, 1/4].
Returns
-------
phi_minus, phi_plus
Two solutions with phi_minus <= K <= phi_plus.
"""
t_minus, t_plus = self.t_solutions(g)
phi_minus = float(self.K) * (t_minus ** (1.0 / float(self.n)))
phi_plus = float(self.K) * (t_plus ** (1.0 / float(self.n)))
return phi_minus, phi_plus
[docs]@dataclass(frozen=True)
class CaptureIslandWindow:
r"""
Compute capture-island (dome-like) concentration intervals from time-scale windows.
The strict gated window is:
.. math::
\tau_{\text{net}} \lesssim \tau(\phi) \lesssim \tau_{\text{obs}}
Using either:
- Rate-averaged :math:`\tau_{\text{eff}}` (S53) => lambda bounds (S60)
- Time-averaged :math:`\langle\tau\rangle` (S69) => lambda bounds (S70)
With:
.. math::
\lambda(\phi) = \lambda_0 \mathcal{G}(\phi)
Therefore:
.. math::
\mathcal{G}(\phi) \in [g_-, g_+], \quad g_{\\pm} = \lambda_{\\pm} / \lambda_0
Because :math:`G(\phi)` is bell-shaped, the set
:math:`\\{\phi: G \in [g_-, g_+]\\}` is typically two disjoint intervals
flanking the peak (an annulus/band around :math:`\phi \approx K`).
Parameters
----------
lambda0 : float
Geometric/multivalency strength parameter :math:`\lambda_0`.
boundaries : DomePhaseBoundaries
Analytic dome boundary solver for :math:`G(\phi) = g`.
See Also
--------
DomePhaseBoundaries : Solver for bell-shaped gating function boundaries.
"""
lambda0: float
boundaries: DomePhaseBoundaries
[docs] def exists(self, g_minus: float) -> bool:
r"""
Check whether an island exists (i.e., threshold below peak).
Parameters
----------
g_minus : float
Lower gating threshold :math:`g_-`.
Returns
-------
exists : bool
True if :math:`g_- \leq 1/4` (the maximum of the gating function).
"""
return 0.0 <= float(g_minus) <= 0.25
[docs] def g_from_lambda(self, lam: float) -> float:
r"""
Convert lambda threshold to g threshold.
Parameters
----------
lam : float
Lambda threshold :math:`\lambda`.
Returns
-------
g : float
Gating threshold :math:`g = \lambda / \lambda_0`.
Raises
------
ValueError
If `lambda0` is not positive.
"""
if float(self.lambda0) <= 0.0:
raise ValueError("lambda0 must be > 0.")
return float(lam) / float(self.lambda0)
[docs] def phi_intervals_for_g_band(self, g_minus: float, g_plus: float) -> list[tuple[float, float]]:
r"""
Return concentration intervals where :math:`G(\phi)` is between :math:`[g_-, g_+]`.
Parameters
----------
g_minus : float
Lower threshold (>=0).
g_plus : float
Upper threshold (>= g_minus).
Returns
-------
intervals : list of tuple
List of (phi_left, phi_right). Usually two disjoint intervals.
Notes
-----
For :math:`0 < g_- < g_+ \leq 1/4`:
- Outer interval from :math:`g_-`: :math:`[\phi_-(g_-), \phi_+(g_-)]`
- Inner interval from :math:`g_+`: :math:`[\phi_-(g_+), \phi_+(g_+)]`
The band is:
:math:`[\phi_-(g_-), \phi_-(g_+)] \cup [\phi_+(g_+), \phi_+(g_-)]`.
"""
g_minus = float(g_minus)
g_plus = float(g_plus)
if g_minus < 0.0 or g_plus < 0.0 or g_plus < g_minus:
raise ValueError("Require 0 <= g_minus <= g_plus.")
if g_plus == 0.0:
return []
if g_plus > 0.25:
raise ValueError("g_plus exceeds peak 1/4; no real boundaries.")
# If g_minus is zero, the "outer" boundary is at 0 and infinity.
if g_minus == 0.0:
phi_in_minus, phi_in_plus = self.boundaries.phi_solutions(g_plus)
return [(0.0, phi_in_minus), (phi_in_plus, float("inf"))]
phi_out_minus, phi_out_plus = self.boundaries.phi_solutions(g_minus)
phi_in_minus, phi_in_plus = self.boundaries.phi_solutions(g_plus)
# Ensure ordering (numerically robust)
a1, a2 = min(phi_out_minus, phi_in_minus), max(phi_out_minus, phi_in_minus)
b1, b2 = min(phi_in_plus, phi_out_plus), max(phi_in_plus, phi_out_plus)
intervals: list[tuple[float, float]] = []
if a2 > a1:
intervals.append((a1, a2))
if b2 > b1:
intervals.append((b1, b2))
return intervals