Worked Analysis Of A Bistable Feedback Switch

From Feedback Equation to a Switch

Positive feedback can amplify an input: greater pathway activity promotes further activity. If that reinforcement saturates while activity is also removed at a rate proportional to its level, the resulting balance can have more than one stable solution. The system may then remain in a low-activity state or a high-activity state at the same input. Which state it occupies depends on its history, so changing the input upward and downward can produce different switching thresholds.

The model below is a compact example of a bistable feedback switch. It is deterministic and describes a single pathway activity variable; it is not a claim that every pathway with positive feedback is bistable. Its equilibria and switching predictions depend on the specified feedback strength, basal production, and assumptions.

Model: Let $x\geq 0$ be dimensionless pathway activity, $u\geq 0$ a constant external input, $b$ basal production, and $\tau>0$ a time constant. The saturating positive-feedback model is $$ \tau\frac{dx}{dt}=u+b+\frac{Vx^2}{K^2+x^2}-x, $$ where $V>0$ is the maximum feedback contribution and $K>0$ is the activity giving half-maximal feedback. For the worked analysis, use $b=0.15$, $V=1.6$, and $K=1$. Activity, input, and these production terms are scaled consistently as dimensionless quantities; $\tau$ sets the time scale.

The model assumes that input and parameters are fixed while an equilibrium is assessed, feedback follows the specified saturating function, and activity loss is first-order. It ignores stochastic fluctuations, delays, and spatial variation. Because $\tau$ is positive, it affects the speed of change but not the equilibrium locations or their stability classification.

Equilibria and Local Stability

An equilibrium is a value $x^*$ at which activity no longer changes, so the right-hand side of the differential equation is zero. For $K=1$, rearranging that condition gives an input value as a function of equilibrium activity:

u=q(x^*)-b,\qquad q(x)=x-\frac{Vx^2}{1+x^2}.

This form is useful because it describes the equilibrium curve in the input–activity plane. At a specified input, intersections with this curve are candidate steady states. Local stability follows from the slope of the differential equation with respect to activity. If that derivative is negative at an equilibrium, a small displacement produces a restoring change; if it is positive, a small displacement grows.

\left.\frac{d}{dx}\left(u+b+\frac{Vx^2}{1+x^2}-x\right)\right|_{x=x^*} =\frac{2Vx^*}{(1+(x^*)^2)^2}-1.

For $V=1.6$, this derivative is negative when $q'(x^*)>0$, and positive when $q'(x^*)<0$, since the derivative of the differential-equation right-hand side is $-q'(x^*)$. Thus equilibria on increasing portions of the curve $u=q(x)-b$ are locally stable, and those on decreasing portions are unstable. This is a local result: it classifies nearby behavior, not the response to every possible perturbation.

Worked Example: Finding the Equilibria at an Intermediate Input

Set $u=0.027$, with $b=0.15$, $V=1.6$, and $K=1$. The equilibrium equation is $$ 0.177+\frac{1.6x^2}{1+x^2}-x=0, $$ because $u+b=0.027+0.15=0.177$. Multiplying by $1+x^2$ and rearranging gives $$ x^3-1.777x^2+x-0.177=0. $$

Numerically solving this cubic gives $x\approx0.3661$, $x\approx0.5863$, and $x\approx0.8246$. These are three candidate equilibrium activities at the same input. For example, substitution of the first rounded value gives $0.177+1.6(0.3661)^2/[1+(0.3661)^2]-0.3661$, which is approximately zero; the other two roots likewise solve the equation to the displayed precision.

To classify them, evaluate the derivative $3.2x/(1+x^2)^2-1$. At $x\approx0.3661$, this is approximately $-0.089$, so the low-activity equilibrium is stable. At $x\approx0.5863$, it is approximately $+0.039$, so the middle equilibrium is unstable. At $x\approx0.8246$, it is approximately $-0.065$, so the high-activity equilibrium is stable. The two stable states are separated by an unstable state.

Fold Points and Switching Thresholds

A fold occurs where the equilibrium curve turns around. At that point, a stable and an unstable equilibrium meet and disappear as the input is changed further. Differentiate $q(x)$; folds satisfy $q'(x)=0$:

q'(x)=1-\frac{3.2x}{(1+x^2)^2}=0, \qquad (1+x^2)^2=3.2x.

The positive solutions are approximately $x=0.4559$ and $x=0.7187$. The corresponding input is $u=q(x)-0.15$. Keeping enough digits in the calculation matters because these values define the predicted switching thresholds.

Worked Example: Calculating Both Fold Inputs

At the first fold, use the displayed rounded activity $x=0.4559$. Then $x^2=0.207844$, and $$ q(0.4559)=0.4559-\frac{1.6(0.207844)}{1+0.207844} \approx0.180574. $$ Subtracting basal production gives $u\approx0.180574-0.15=0.030574$. To five significant figures, the upper input fold is $u_{\mathrm{up}}\approx0.03057$.

At the second fold, $x=0.7187$ gives $x^2=0.51652969$, so $$ q(0.7187)=0.7187-\frac{1.6(0.51652969)}{1+0.51652969} \approx0.173740. $$ Thus $u\approx0.173740-0.15=0.023740$, or $u_{\mathrm{down}}\approx0.02374$ to five significant figures. The small difference between the two thresholds is real: the upper fold is about $0.03057$, and the lower fold is about $0.02374$.

The activity values used here are rounded approximations to the fold locations, so the computed inputs are approximate as well. Reporting $q$ and $u$ to six decimal places is consistent with evaluating the displayed rounded activities; the thresholds are best interpreted at the precision supported by those activities.

Between the fold inputs, the equilibrium curve has three intersections: low and high stable states with an unstable state between them. Below the lower fold, only the low-activity stable branch remains; above the upper fold, only the high-activity stable branch remains. At a fold itself, the two meeting equilibria are marginal in the local linearization, so it is not appropriate to describe the fold equilibrium as an ordinary, robust stable state.

The labels $u_{\mathrm{up}}$ and $u_{\mathrm{down}}$ refer to the direction of an input ramp, not to the activity at the fold. During an upward ramp, the system can follow the low stable branch until it reaches the upper input fold. During a downward ramp from the high state, it can follow the high branch until it reaches the lower input fold. These distinct thresholds are the hysteresis predicted by this deterministic model.

Using the Equilibrium Curve to Predict a Ramp

A numerical simulation can show the time course as well as the equilibria. For a fixed input and a time step $\Delta t$, a simple forward Euler update is

x_{j+1}=x_j+\frac{\Delta t}{\tau}\left(u+b+\frac{Vx_j^2}{1+x_j^2}-x_j\right).

This approximation requires a sufficiently small time step for the chosen $\tau$ and dynamics. To simulate a changing input, update $u$ according to a specified ramp and apply the activity update at each time step. The result should be checked by reducing $\Delta t$; a qualitative switch that vanishes or shifts substantially with a smaller step may be a numerical artifact.

Worked Example: Predicting an Upward and Downward Ramp

Suppose the model is initialized on its low-activity branch at $u=0.01$, and the input is then raised slowly to $0.04$ before being lowered slowly back to $0.01$. These inputs span both calculated folds: $0.01<0.02374<0.03057<0.04$.

On the upward ramp, the system follows the low stable branch while it exists. Near $u_{\mathrm{up}}\approx0.03057$, that branch meets the middle unstable branch at activity $x\approx0.4559$. Above the fold, those equilibria no longer exist, so the trajectory moves toward the remaining high-activity stable state. At the fold input, the cubic equilibrium equation has a repeated root near $0.4559$; its third root is near $0.8688$, the high-activity equilibrium available after the low branch disappears.

On the downward ramp, the activity follows the high branch instead. Near $u_{\mathrm{down}}\approx0.02374$, the high stable branch meets the middle unstable branch at $x\approx0.7187$. The low-activity equilibrium at this fold is near $0.3363$, the third root of the cubic at that input. Once the high branch disappears, activity moves toward the low branch. Thus, at an input such as $u=0.027$, the predicted state depends on the direction of approach: it can be near $0.3661$ after an upward-starting low-state history or near $0.8246$ after a downward-starting high-state history.

A finite ramp need not switch at exactly the quasi-static fold values. The trajectory takes time to respond, and a ramp that is not sufficiently slow can produce delay or lag. Therefore these thresholds are predictions for the equilibrium structure; simulated or observed switching during a finite-rate ramp should be interpreted in light of the ramp speed and time scale.

Checking the Analysis and Its Limits

Equilibrium calculations, local stability analysis, and dynamic simulation answer related but different questions. Solving the steady-state equation identifies candidate states. The derivative classifies their local stability. Simulation shows a trajectory under a chosen initial condition and input schedule. Agreement among the three is a useful check, but none replaces the others.

Check the root count. In the bistable input interval, a numerical solver should find three nonnegative equilibrium roots, not just one. A single root may mean the solver used an unsuitable starting point or that only one branch was searched.
Check the derivative sign. The outer equilibria should have negative differential-equation derivative and the middle one a positive derivative. If these signs are reversed, recheck the derivative or the stability convention.
Check the ramp history. At a fixed input in the bistable interval, distinct final states require initial conditions or ramp histories that place the system in different basins of attraction. One simulation from one initial state cannot establish that both stable states exist.
Check numerical sensitivity. Repeat simulations with smaller time steps and slower ramps. Finite-step errors and finite-rate lag can shift apparent transition points relative to equilibrium folds.

The equilibrium curve is often clearer than a long simulation when the question is how many steady states exist and where folds occur. Simulation is needed to examine transients, ramp-rate effects, and trajectories from specified initial conditions. Alternative parameter estimation or more detailed pathway models may be warranted when data show delays, stochastic switching, or feedback kinetics inconsistent with the simple saturating term; those behaviors are not captured by this one-variable deterministic model.

Common Pitfall

Common mistake: “The response is steep, so it must be bistable.” A steep but single-valued equilibrium curve can still have only one stable state at each input. Bistability requires two stable equilibria at the same input, separated here by an unstable equilibrium; hysteresis requires different switching thresholds on upward and downward input changes. Verify the number of roots, their stability, and the fold locations rather than inferring bistability from a steep response alone.

The model also does not imply that a pathway will remain in a state indefinitely in a noisy biological system. Noise can cause transitions between states, and finite observation windows can affect whether switching is seen. The deterministic fold thresholds describe where equilibria appear or disappear in the stated equations, not guaranteed thresholds for every cell or experiment.

Key takeaway: For this saturating feedback model, the equilibrium curve has two stable branches over an input interval bounded approximately by $u_{\mathrm{down}}=0.02374$ and $u_{\mathrm{up}}=0.03057$. A slow upward ramp is predicted to switch from low to high activity at the upper fold; a slow downward ramp is predicted to switch from high to low activity at the lower fold.

Check Your Understanding

Try each question, then compare your reasoning with the solutions below.

At $u=0.027$, how many equilibria does the model have, and which are locally stable?
Why does the upward switching threshold differ from the downward switching threshold in this model?
At a fixed input within the bistable interval, what must be specified to predict which stable state a simulation reaches?
What does reducing the time step check in a numerical ramp simulation, and what does it not eliminate?

Solutions. (1) The cubic has three nonnegative roots, approximately $0.3661$, $0.5863$, and $0.8246$. The outer roots are locally stable because the differential-equation derivative is negative there; the middle root is unstable because it is positive. (2) The low and high stable branches terminate at different folds: the low branch ends at the upper input fold, while the high branch ends at the lower input fold. (3) The initial activity or the input history must be specified. In the bistable interval, one stable equilibrium does not uniquely determine the final state. (4) Reducing the time step checks sensitivity to discretization error in the numerical integration. It does not remove finite-rate lag, stochastic effects absent from this deterministic model, or the model’s structural assumptions.