Saddle-Node Bifurcations And Hysteresis In Pathways

From bistability to switching thresholds

Bistability means that, over some range of conditions, a pathway can settle into either of two stable activity states. It does not by itself specify when the pathway switches as conditions change. That requires locating the saddle-node bifurcations at which a stable equilibrium and an unstable equilibrium meet and disappear. When the input is ramped slowly, the two folds can define different thresholds for switching upward and downward. The resulting dependence on the direction of the input change is hysteresis.

We use a dimensionless autoactivation model with pathway activity \(x\geq 0\) and external input \(u\):

\frac{dx}{dt}=f(x;u)=u+\frac{A x^2}{1+x^2}-x.

Here \(A>0\) is the strength of positive feedback, \(u\) is a constant input during each equilibrium calculation, and the linear loss rate has been scaled to one. Thus time is in units of the loss timescale, while \(x\), \(u\), and \(A\) are dimensionless in this normalized model. The saturating term is a Hill function with exponent two. We treat \(A\) as fixed while varying \(u\).

Definition: An equilibrium \(x^*\) satisfies \(f(x^*;u)=0\). In this one-dimensional model, it is locally asymptotically stable if \(f_x(x^*;u)<0\), and unstable if \(f_x(x^*;u)>0\). A saddle-node (fold) occurs when an equilibrium satisfies both \(f(x;u)=0\) and \(f_x(x;u)=0\), with the equilibrium branch turning as the parameter changes.

This is a deterministic, well-mixed model with fixed parameters and no random fluctuations. Interpreting the folds as observed switching thresholds additionally assumes that input changes slowly enough for the pathway to track a stable equilibrium, and that the system remains near that equilibrium until it loses stability. Rapid ramps, noise, unmodeled feedback, or substantial parameter drift can change observed switching.

Deriving the fold curve

Differentiate the rate equation with respect to activity:

f_x(x;u)=\frac{2Ax}{(1+x^2)^2}-1.

The fold condition \(f_x=0\) gives the feedback strength required for a fold at activity \(x>0\):

A_{\mathrm{fold}}(x)=\frac{(1+x^2)^2}{2x}.

At a fold, substitute this expression into the equilibrium equation. Since \(A_{\mathrm{fold}}x^2/(1+x^2)=x(1+x^2)/2\), the input at the fold is

u_{\mathrm{fold}}(x) =x-\frac{x(1+x^2)}{2} =\frac{x(1-x^2)}{2}.

This substitution matters: the feedback term at the fold is not simply \(x/2\). The fold activity \(x\) must satisfy the feedback-strength equation for the chosen \(A\); the corresponding \(u_{\mathrm{fold}}\) then gives its input threshold.

A short derivative calculation shows when two fold activities exist. Rewrite the fold curve as \(A_{\mathrm{fold}}(x)=1/(2x)+x+x^3/2\). Then

\frac{dA_{\mathrm{fold}}}{dx} =-\frac{1}{2x^2}+1+\frac{3x^2}{2}=0 \quad\Longleftrightarrow\quad 3x^4+2x^2-1=0.

The positive solution is \(x=1/\sqrt{3}\), and the fold curve has its minimum there: \(A_{\mathrm{crit}}=8\sqrt{3}/9\approx 1.5396\). For \(A>A_{\mathrm{crit}}\), there are two positive fold activities. For both folds to occur at positive inputs in this model, the upper fold activity must be below one; because \(A_{\mathrm{fold}}(1)=2\), choosing \(A_{\mathrm{crit}}<A<2\) gives two positive input thresholds. This range is useful for illustrating a hysteresis loop under nonnegative input.

Worked Example: Evaluating a Fold-Curve Point

Evaluate the fold-strength curve at activity \(x=0.2091\), and distinguish that rounded evaluation from one using a more precise activity.

Substitute the displayed activity into the formula:

A_{\mathrm{fold}}(0.2091) =\frac{(1+0.2091^2)^2}{2(0.2091)} =\frac{(1.04372281)^2}{0.4182} \approx\frac{1.089357}{0.4182} \approx 2.60487 \approx 2.605.

Thus \(2.605\), not \(2.604\), is the result rounded to three decimal places for the displayed substitution. If a fold activity is instead given more precisely as approximately \(0.2091488\), evaluating the same curve gives approximately \(2.6044\), which is about \(2.604\) to three decimal places. A reported threshold must be consistent with the precision of the activity actually substituted.

Locating both thresholds in a positive-input loop

For a concrete hysteresis curve, set \(A=A_{\mathrm{fold}}(0.35)\). This selects a feedback strength just above the critical value and makes one fold activity exactly \(0.35\). The second fold is found by solving \(A_{\mathrm{fold}}(x)=A\). Both folds are calculated below; the lower-activity fold gives the larger input threshold in this example.

Worked Example: Calculating the Two Fold Thresholds

First calculate the fixed feedback strength:

A=\frac{(1+0.35^2)^2}{2(0.35)} =\frac{(1.1225)^2}{0.7} =\frac{1.26000625}{0.7} \approx 1.8000089.

The fold equation can be written as \(x^4+2x^2+1-2Ax=0\). With this \(A\), its two positive roots are approximately \(x_{\mathrm{low}}=0.3500\) and \(x_{\mathrm{high}}=0.88725\). Substitution into the fold-input formula gives

u_{\mathrm{fold}}(0.3500) =\frac{0.3500(1-0.3500^2)}{2} =\frac{0.3500(0.8775)}{2} =0.1535625,

u_{\mathrm{fold}}(0.88725) =\frac{0.88725(1-0.88725^2)}{2} \approx \frac{0.88725(0.21280)}{2} \approx 0.09440.

The thresholds are therefore approximately \(u_{\mathrm{down}}=0.09440\) and \(u_{\mathrm{up}}=0.1536\). At the upper threshold, the low-activity stable equilibrium meets the intermediate unstable equilibrium; increasing input past this fold removes the low state and produces an upward switch. At the lower threshold, the high-activity stable equilibrium meets the intermediate unstable equilibrium; decreasing input past this fold removes the high state and produces a downward switch. The two thresholds are distinct, and both are positive.

Stability and the hysteresis curve

The derivative \(f_x\) determines the local stability of each equilibrium. At very small positive \(x\), \(f_x\) is negative. Between the two fold activities, it is positive, and above the second fold activity it is negative again. Consequently, when three equilibria exist, the low- and high-activity equilibria are stable and the middle equilibrium is unstable. At either fold, \(f_x=0\), so linear stability alone does not classify the equilibrium at the threshold; the branch geometry identifies the collision.

The equilibrium curve can also be obtained by solving \(f=0\) for the input:

u_{\mathrm{eq}}(x)=x-\frac{Ax^2}{1+x^2}.

Plotting this curve against \(x\) gives a turning, or S-shaped, equilibrium diagram when \(A>A_{\mathrm{crit}}\). Its turning points are the folds. To construct the hysteresis path, follow the stable lower branch as \(u\) rises, switch to the stable upper branch at \(u_{\mathrm{up}}\), and then follow that branch as \(u\) falls. The downward switch occurs at \(u_{\mathrm{down}}\). The unstable middle branch is part of the equilibrium diagram, but a deterministic system does not remain on it under small perturbations.

Fix the feedback strength.
Choose \(A\) and verify that it exceeds \(A_{\mathrm{crit}}\). Check whether both fold inputs are in the input range being studied.

Find fold activities and inputs.
Solve \(A=A_{\mathrm{fold}}(x)\) for positive fold activities, then calculate \(u_{\mathrm{fold}}=x(1-x^2)/2\) at each root.

Classify equilibria and trace the path.
Use the sign of \(f_x\) on each branch. Follow a stable branch as input changes and switch when that branch terminates at a fold.

Worked Example: Tracing an Upward and Downward Ramp

Use the thresholds from the preceding example. Suppose input is increased slowly from \(u=0.08\) to \(u=0.17\), then decreased slowly to \(u=0.08\). Assume the system begins on the low-activity stable branch and tracks stable equilibria adequately between thresholds.

During the upward ramp, the pathway remains on the low branch through \(u=0.09440\). Between \(0.09440\) and \(0.1536\), both stable states exist, but the pathway’s current state determines which one it follows. At \(u=0.1536\), the low branch ends in a fold; the pathway switches to the high branch. The ramp reaches \(u=0.17\) on that high branch.

As input is reduced, the pathway remains on the high branch through \(u=0.1536\), even though the low state is also available. It switches downward only when the high branch ends at \(u=0.09440\). Thus the activity at an input between the thresholds depends on the direction of the ramp and the state occupied previously. This path dependence is the hysteresis loop; it is not a claim that the equilibrium at a given input is uniquely determined by input alone.

Diagnostics, alternatives, and limitations

A useful diagnostic is to calculate equilibria at several inputs on both sides of each proposed fold and check that the number and stability of the roots change as expected. A numerical equilibrium solver should use multiple initial guesses: a single initial condition can find one stable root while missing another stable root and the unstable root. Near a fold, roots approach one another, so coarse input grids and loose numerical tolerances can obscure the turning point. Continue the equilibrium curve in activity or use a fold-solving method based on both \(f=0\) and \(f_x=0\), rather than identifying thresholds only from a visibly abrupt jump in a time simulation.

The equilibrium diagram and a time-course simulation answer related but different questions. The diagram locates all equilibria, including unstable ones, and exposes both folds. A simulation shows the trajectory under a specified input schedule and initial condition, but it can miss branches it never visits. Directly sweeping input in simulations is informative as a check, not a substitute for locating equilibria. A rapid ramp can cause delayed switching beyond a static fold; the quasi-static thresholds derived here are equilibrium thresholds, not universal measured thresholds.

The Hill model is a compact representation, not a detailed biochemical mechanism. Its thresholds depend on the assumed feedback function, normalized loss, and fixed feedback strength. If the chosen \(A\) yields a fold at a negative input, that fold is mathematically present but is not a positive-input switching threshold for an experiment restricted to \(u\geq0\). Also, deterministic hysteresis should not be confused with a noise-induced transition: random switching can occur before a deterministic branch disappears, a question requiring a stochastic model.

Using only equilibrium roots at one input. A single calculation cannot establish two switching thresholds; locate both folds and check their input values.
Calling every equilibrium stable. Evaluate \(f_x\). The intermediate branch in the three-equilibrium region is unstable.
Swapping the thresholds. In the positive-input example, the lower-activity fold has the higher input threshold. Determine the direction of switching from which stable branch terminates, not from fold activity alone.
Confusing a static fold with a measured jump. Finite ramp speed, noise, and observation resolution can shift or blur an observed transition.

Key takeaway: For the Hill autoactivation model, solve \(A=A_{\mathrm{fold}}(x)\) to locate fold activities, then calculate \(u_{\mathrm{fold}}=x(1-x^2)/2\). The stable branches and their distinct fold inputs determine the upward and downward paths of the hysteresis curve.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions.

For \(A=1.8\), the model has two positive fold activities. What equations must be satisfied at each fold?
Why is \(f_x=0\) alone insufficient to identify a fold of the equilibrium curve?
In the worked ramp, at an input between \(0.09440\) and \(0.1536\), why can two different pathway activities be observed?
A simulation finds one stable equilibrium at a chosen input. Does that rule out bistability? Explain.

Solutions. (1) Each fold satisfies \(f(x;u)=0\) and \(f_x(x;u)=0\); equivalently, \(A=A_{\mathrm{fold}}(x)\) and \(u=x(1-x^2)/2\). (2) \(f_x=0\) identifies a stationary point of the rate as a function of activity, but without \(f=0\) that point is not an equilibrium. (3) The pathway can occupy either stable branch depending on its prior state and the direction of the input ramp; the intermediate equilibrium is unstable. (4) No. A single initial condition or solver start may converge to just one stable root. Multiple initial guesses and an equilibrium or fold analysis are needed to check for other roots.