From model structure to a defensible behavior claim

A worked method for deciding what a reaction-model family permits, what one parameter choice produces, and which calculation supports each conclusion.

CCBS · EST 5053 · Fall 2026 · lecture 4 · lecturer-written exposition · study alone after the core lecture

Lecturer-written companion, released 15 September 2026 after the September 10 lecture. Álvaro’s student slides are linked separately on the course schedule. The lineage traces the methods, and the audience extension develops the Week 2 student contributions.

The core lecture asks for a map from biochemical structure to dynamical function. This exposition teaches how to build one without confusing a displayed trajectory with a conclusion about a design. A reaction mechanism, admissible state space, kinetic assumptions and parameter ranges define a model family. A numerical run chooses one member of that family and one initial state. Global geometry, local linearization and scalar certificates can instead establish statements that survive specified changes within the family.

1Choose a state

State is not merely what can be measured. It is the present information required by the model to determine the instantaneous future.

Consider a mass attached to a spring. If only position $q$ is recorded, the same value of $q$ can occur while the mass moves left or right. Position therefore does not determine $\dot q$. Velocity must be included. Starting from

Here $m>0$ is mass, $k>0$ is spring stiffness and $c\ge0$ is the damping coefficient. Define $v=\dot q$ and the state vector $z=(q,v)$. Then

This is an exact change of representation, not an approximation. The second-order scalar equation and the two-dimensional first-order equation contain the same information. The benefit is conceptual: equation (1) assigns one vector to every point in the $(q,v)$ plane. Trajectories cannot cross there, because crossing would assign two future velocities to one state.¹

The same state test applies to cells. Protein concentration alone may not determine its future if a promoter can be in two configurations with different production rates. One can include promoter state explicitly, include an external input as $F(x,t)$, or derive an effective law by eliminating a fast hidden state. The choice changes what the model claims.

Choose variables that determine the model's instantaneous future, and distinguish state augmentation from an approximation.

2Convert memory into state

Higher order, delay, and hidden intermediates all express memory. State augmentation makes that memory explicit.

An $n$th-order scalar equation can be written as $n$ first-order equations. For

set $z_1=q,z_2=\dot q,\ldots,z_n=q^{(n-1)}$. The resulting companion system is

A finite chain of maturation steps works the same way. If immature protein $P_0$ matures through $P_1$ to fluorescent $P_2$, then the current value of $P_2$ does not reveal the pipeline already committed upstream. Retaining $(P_0,P_1,P_2)$ turns the apparent memory into ordinary first-order flow.

A true fixed delay is different: its exact state is a history function, hence infinite-dimensional. A chain of intermediate states can approximate a distributed delay, but the approximation and its timescale must be declared. “Everything is first order” means that an appropriate state representation exists, not that one measured variable always suffices.

Rewrite a finite-order equation as a first-order state model, and recognize when a delay requires a history rather than finitely many variables.

3Prove the phase-line rules

For one autonomous variable, ordering makes the sign of one function almost a complete global analysis.

Let $\dot x=f(x)$, with $f$ locally Lipschitz so solutions are locally unique. A fixed point satisfies $f(x^*)=0$. Between two consecutive roots, continuity prevents the sign of $f$ from changing. Therefore every nonstationary solution in that interval is monotone.

If $f$ is twice continuously differentiable near a simple root, Taylor expansion gives $f(x^*+u)=f'(x^*)u+O(u^2)$. To first order, $\dot u=f'(x^*)u$, so $u(t)\approx u(0)e^{f'(x^*)t}$. A negative derivative attracts from both sides, while a positive derivative repels. When $f'(x^*)=0$, inspect higher-order terms or the exact signs of the field.

The equation $\dot x=x^2$ shows why the three views are not redundant. Its origin attracts from the left and repels from the right. Separation gives

The phase line predicts the direction, but only equation (3) shows that a positive initial condition escapes in finite time $1/x_0$. Conversely, ordering proves a result that no finite collection of traces can prove: a scalar autonomous system has no nonconstant periodic orbit. A periodic trajectory would have to return to a previous value after moving monotonically, contradicting uniqueness.

Use phase-line signs to establish direction and basins, then solve or estimate the travel time when finite-time escape matters.

4Follow a saddle-node family

A bifurcation is a change in the organization of a family of vector fields, not simply a curve drawn in parameter space.

Consider the saddle-node normal form.²

Solving $f=0$ gives no real fixed point for $r<0$, one double root at $r=0$, and two roots $x^*_{\pm}=\pm\sqrt r$ for $r>0$. Since $f_x=-2x$, the positive branch is stable and the negative branch unstable. At the collision, $f_x=0$, so the hyperbolic linear test fails exactly where the topology changes.

The local return rate also becomes measurable. Write $x=\sqrt r+u$ near the stable branch:

Thus $\tau_{\mathrm{relax}}=(2\sqrt r)^{-1}$ diverges as the collision is approached. This is the mechanism behind critical slowing near a generic fold. In data, increased recovery time can be confounded by noise color, nonstationary forcing, sparse sampling, or a transition that is not a fold, so it is evidence to test rather than a universal alarm.⁴

Distinguish an equilibrium branch from a trajectory and derive the return-time divergence near a generic fold.

5Open the toggle circuit

Mutual repression is a network of regulated addition and removal, not an arbitrary pair of sigmoids.

Let gene $G_X$ produce repressor $X$ and gene $G_Y$ produce repressor $Y$. Each protein represses the other gene. At a composite level:

The removal processes need not have the same molecular origin. If $N_X$ is the molecule count and a bacterial cell grows as $V(t)=V_0e^{\lambda t}$, a protein that is stable as a molecule can still be diluted as a concentration. With production $dN_X/dt=\rho_XV$ and concentration $X=N_X/V$,

Molecular degradation adds its own first-order contribution, so an effective removal rate can be $\delta_X=\lambda+k_{\mathrm{deg},X}$. This connects the circuit directly to reaction and growth bookkeeping.

Separate regulated production, molecular degradation and growth dilution before writing a gene circuit's concentration equations.

6Derive the effective repression law

A Hill term is a reduced occupancy law. Its assumptions are part of the model and should travel with the formula.

For the $X$-producing promoter, let $G_X$ denote the active configuration and $C_{G_XY}$ a configuration containing that promoter and repressor $Y$. Suppose fast binding can be summarized by

Eliminate $C_{G_XY}$ using promoter conservation:

If promoter redistribution is fast relative to protein accumulation and production is proportional to active occupancy, the slow protein equations become

Equation (8) is an effective binding relation. The exponent can encode concerted switching, cooperative occupancy, or a phenomenological fit. It is not automatically the molecularity of one elementary collision. Independent identical sites with production only when all are empty would instead give $(1+Y/K)^{-n}$. Sequential binding, unequal sites, sequestration, and nonequilibrium promoter cycles yield still other rational functions. The reduction is useful because it states what is being compressed.

Derive an effective repression law from an explicit occupancy assumption and name the conditions under which the eliminated states can be ignored.

7Nondimensionalize and find every state

Scaling separates circuit shape from units and makes a complete equilibrium calculation possible.

Let $x=X/K_X$, $y=Y/K_Y$, and use a common removal time $\tau=\delta t$ for a symmetric teaching case. Define $a_X=\alpha_X/(\delta K_X)$ and $a_Y=\alpha_Y/(\delta K_Y)$. With $a_X=a_Y=4$ and $n_X=n_Y=2$,

At a symmetric equilibrium, $x=y=s$ and $s^3+s-4=0$, whose positive root is $s\approx1.3788$. The asymmetric solutions are exact:

or approximately $(0.2679,3.7321)$ and $(3.7321,0.2679)$. Substitution into both equations is essential. A point on one nullcline stops only one component and is generally still moving.

Nondimensionalize a circuit and verify every proposed equilibrium in both component equations.

8Assemble basins from nullclines

Nullclines locate the component sign changes. Their intersections and the flow between them organize the decision.

The $x$-nullcline is $x=4/(1+y^2)$. Below it, $f>0$ and motion has a rightward component; above it, $f<0$. The $y$-nullcline is its reflection. Testing one point in each region determines the arrow directions without integrating.

Symmetry gives a useful proof. If $x=y$, then $f(x,x)=g(x,x)$, so $d(x-y)/dt=0$. The diagonal is invariant. Uniqueness prevents a trajectory from crossing it, and each side flows to its corresponding asymmetric expression state. Thus the diagonal is the stable manifold of the middle saddle and the basin boundary.

A steep input-output curve alone is not bistability. Bistability requires two stable invariant outcomes at the same parameter value, separated by a threshold. Hysteresis additionally concerns how the state follows branches when a parameter is swept. Gardner, Cantor, and Collins constructed a mutually repressive switch in Escherichia coli; the dimensionless system here is a transparent teaching model of the same dynamical design, not a parameter fit to that construct.³

Use nullcline signs and invariant sets to distinguish a basin boundary from a single zero-velocity component.

9Derive the linearized system

The Jacobian matrix becomes useful only after it is evaluated at a fixed point and placed inside a displacement equation.

For $\dot x=f(x,y)$, $\dot y=g(x,y)$, let $(x^*,y^*)$ be a fixed point and define the displacement vector $u=(u_1,u_2)=(x-x^*,y-y^*)$. Taylor expansion gives

The canonical linearized system is the separate first-order model

For the toggle,

Every sign has biological content. The diagonal entries are effective protein removal. Growth dilution alone contributes $-\lambda_i$ to the corresponding diagonal before other interactions are added. The off-diagonal entries are negative because each protein represses production of the other. Evaluating at different intersections gives different linear systems, so “the Jacobian of the model” is incomplete language unless the state is named.

Define displacement from a specified equilibrium and write the Jacobian inside the linearized dynamical system.

10Read modes and the Hurwitz rule

Eigenvalues are rates because an eigendirection turns a vector equation into one scalar exponential.

If $Aw=\lambda w$ and $z(t)=a(t)w$, then

For a real eigenvalue, the sign controls growth or decay along that eigendirection. For $\lambda=\sigma+i\omega$, the real part $\sigma$ controls the envelope $e^{\sigma t}$ and $\omega$ controls rotation. A general solution is a sum of modes, with polynomial factors added when a matrix has nontrivial Jordan blocks.

The intuition is eventual decay. Every possible initial displacement is assembled from modal components. If even one component grows, a small perturbation with a projection onto it eventually leaves. If all exponential envelopes decay, polynomial prefactors cannot defeat them and every linear trajectory tends to the equilibrium. The distance need not decrease at every instant. For the nonlinear system, the remainder is smaller than linear near the equilibrium, so the strict linear stability result persists in a sufficiently small neighborhood.

Use eigenvalue real parts to test local attraction, and keep the nonlinear scope distinct from the global linear result.

11Use trace and determinant

For a real two-dimensional system, two scalar combinations classify all hyperbolic linear portraits.

Since $\lambda_++\lambda_-=\tau$ and $\lambda_+\lambda_-=\Delta$, the two-dimensional Hurwitz condition is

A negative determinant forces two real eigenvalues of opposite sign, hence a saddle. With positive determinant, negative trace puts both real roots on the negative side or gives a conjugate pair with negative real part. The discriminant $D=\tau^2-4\Delta$ distinguishes nodes from spirals but does not decide stability by itself.

For the two asymmetric toggle states, equation (15) has eigenvalues $-1/2$ and $-3/2$, so both are stable nodes. At the symmetric state $s\approx1.3788$, each off-diagonal entry is $-s^3/2\approx-1.3106$. The eigenvalues are therefore $-1\pm1.3106$, one positive and one negative. This is the saddle whose stable direction is tangent to the basin boundary.

The tangent is not the full separatrix. Nor does a stable Jacobian prove that every initial condition reaches that equilibrium. Those are global questions, which require invariant regions, nullclines, monotonicity, a Lyapunov argument, or another global tool.

Apply trace and determinant to classify planar linearizations without mistaking their local conclusions for a global basin argument.

12Know when first order is silent

Purely imaginary eigenvalues report the absence of first-order radial motion, not the presence of a nonlinear center.

Consider three systems written in polar coordinates:

At the origin, all three have the same Cartesian Jacobian $\left[\begin{smallmatrix}0&-1\\1&0\end{smallmatrix}\right]$ and eigenvalues $\pm i$. Yet the higher-order radial laws make trajectories spiral inward, stay on circles, or spiral outward.

This is why the Hurwitz theorem uses strict inequalities. If an eigenvalue lies on the imaginary axis, the nonlinear remainder is no longer dominated by first-order contraction. One must calculate the relevant nonlinear terms, find a conserved quantity, construct a Lyapunov function, or apply a bifurcation theorem.

Recognize when zero real parts leave the nonlinear stability question open and choose a higher-order or global test.

13Make energy a certificate

A Lyapunov function is a generalized energy: positive away from the target and nonincreasing along the actual vector field.

For the spring state in equation (1), define

Differentiate along trajectories:

When $c=0$, equation (21) gives $\dot E=0$. Every nonzero start remains on its energy ellipse, so the origin is stable but not attracting. The purely imaginary eigenvalues are consistent with this conclusion, but energy proves it globally for the linear spring.

When $c>0$, $E$ is positive definite and nonincreasing. The derivative is not strictly negative: it vanishes on the entire line $v=0$. To finish correctly, inspect that set. If $v=0$ but $q\ne0$, then $\dot v=-(k/m)q\ne0$, so the trajectory immediately leaves. The largest invariant set inside $\{\dot E=0\}$ is only $(0,0)$. LaSalle's invariance principle then gives convergence to the origin.

Positive definiteness and decreasing derivative are local unless more is shown. A global conclusion generally needs a forward-invariant domain and compact sublevel sets, often supplied by radial unboundedness. A Lyapunov function proves a statement about trajectories. Merely calling a scalar “energy-like” proves nothing. Modern biochemical-network methods construct such certificates from network structure rather than guessing a mechanical energy.⁷⁸

Check the invariant zero-derivative set and confinement conditions before inferring convergence from an energy balance.

14Work the nonlinear pendulum

The pendulum shows what local linearization captures, what global energy adds, and why the angle's topology matters.

Measure angle $\theta$ from the downward vertical and scale time so $g/L=1$. The ideal pendulum is

Its energy per unit $mL^2$ is

At $(\theta,\omega)=(2\pi n,0)$, the Jacobian is $\left[\begin{smallmatrix}0&1\\-1&0\end{smallmatrix}\right]$, so small oscillations have eigenvalues $\pm i$ and period $2\pi$. Energy supplies the missing nonlinear result: nearby sublevel sets surround the equilibrium, so it is Lyapunov stable, while conservation prevents attraction.

At $((2n+1)\pi,0)$, the lower-left Jacobian entry changes to $+1$, giving eigenvalues $\pm1$ and a saddle. The energy there is $H=2$. Levels below 2 librate around a downward equilibrium, the level 2 contains the separatrices, and levels above 2 rotate. The small-angle approximation captures the local period and ellipses, but it cannot display the upper saddle, the separatrix, amplitude-dependent period, or rotation.

Because $\theta$ and $\theta+2\pi$ are the same physical orientation, the natural phase space is a cylinder rather than an infinite plane. Repeated copies in the figure are a drawing convention, not distinct pendulums.

Use pendulum energy and the cylindrical state space to distinguish oscillation, rotation and the separatrix beyond the small-angle approximation.

15Stabilize the inverted pendulum

The same state-space and Hurwitz tools turn an unstable mechanical equilibrium into a transparent feedback design.

Let $\phi$ be angle from the upright position, and let $a$ be a signed horizontal pivot acceleration chosen so that

Near $\phi=0$, with $\omega=\dot\phi$,

Without control the eigenvalues are $\pm\sqrt{g/L}$, so the upright equilibrium is a saddle. Passive velocity damping alone changes the trace but leaves a negative determinant, hence it cannot remove the unstable direction. Use local proportional-derivative feedback

The closed-loop linearized matrix is

Hurwitz therefore requires $k_d>0$ and $k_p>g$. Proportional action stronger than gravity flips the effective stiffness; derivative action supplies contraction. With $k_d=0$ the linearized system is a center, not an attractor. For the smooth nonlinear model, strict Hurwitz stability proves local exponential stability. It does not prove global recovery because the trigonometric dynamics, finite track, actuator saturation, and unmodeled cart dynamics limit the basin.

Design a locally stabilizing feedback from trace and determinant, then state which nonlinear and actuator limits remain outside that guarantee.

16Separate a center from a limit cycle

Closed trajectories can mean a continuum of conserved levels or one isolated attracting orbit. The distinction is dynamical, not visual.

The ideal spring and pendulum have centers. Every nearby energy value produces a different closed orbit. Perturbing one orbit moves the trajectory onto another; it does not return to the first. By contrast,

has an isolated periodic orbit at $r=1$. Starts with $0<r<1$ move outward and starts with $r>1$ move inward. Both converge to the same orbit, making it a stable limit cycle.

For a smooth planar flow, the Poincaré-Bendixson theorem gives a proof pattern. Find a compact forward-invariant region, exclude fixed points from the limiting set, and conclude that a periodic orbit remains. It does not apply directly to a three-state repressilator or to a delay equation. A closed-looking numerical trace alone also does not invoke the theorem. Vary the step size and initial conditions, integrate long enough to remove transients, and identify the trapping and exclusion arguments.

Distinguish an isolated attracting cycle from a conserved family and state the hypotheses needed for a planar recurrence theorem.

17Find the Lotka-Volterra integral

The classic predator-prey model looks oscillatory, but its hidden conserved function shows that the orbits form a center rather than an attracting clock.

The positive equilibrium is $(x^*,y^*)=(\gamma/\delta,\alpha/\beta)$. Define

Since $u-u^*-u^*\ln(u/u^*)\ge0$ for $u>0$, $V$ is positive relative to the equilibrium. Its derivative cancels exactly:

Thus the positive equilibrium is surrounded by a family of conserved closed curves. The Jacobian also has purely imaginary eigenvalues, but equation (30) supplies the nonlinear reason. The model is not a robust oscillator: changing the initial condition changes the orbit and period, and common ecological modifications destroy the integral.

Differentiate a candidate first integral to test whether a biological oscillation is neutral rather than attracting.

18See how saturation creates a cycle

A saturating predation rate changes the fixed point from a center to an unstable spiral and permits an isolated attracting orbit.

Replace linear prey growth with logistic growth and let predation saturate. Consider the following illustrative system.

Use $\alpha=\gamma=1$, $\beta=\delta=3$, and $K=3$. The positive fixed point is

Here $\tau=1/9>0$ and $\Delta=5/9>0$, with eigenvalues $1/18\pm i\sqrt{179}/18$. The coexistence point is an unstable spiral. The positive quadrant is invariant. Logistic prey growth bounds $x$, and $d(x+y)/dt=2x-x^2/3-(x+y)$ bounds the total. The computed trajectories in Figure 11 approach the same orbit from both sides. A Poincaré-Bendixson proof would additionally need an interior trapping argument excluding equilibria from the limiting set. That final global proof is not supplied by the local eigenvalues or by the figure alone.

Separate a computed attracting cycle from a local instability and from a complete global existence or uniqueness proof.

19Solve two reaction-network portraits

Small mass-action networks connect bifurcation and basin geometry directly to stoichiometric events.

A. A square-root response is a saddle-node slice

Treat $X$ as a clamped input and let the two reactions be

The rate constants are independent and positive. Choose each as one in this scaled illustration. The resulting scaled rates $x$ and $y^2$ give

For $x>0$, the physical equilibrium is $y^*=\sqrt x$, with derivative $-2\sqrt x<0$. The other algebraic root lies outside the concentration domain. As $x\downarrow0$, the return time diverges as $1/(2\sqrt x)$. This is the saddle-node normal form restricted by nonnegativity, and it shows why the state domain belongs in the analysis.

B. A mass-action winner-take-all network

For the chosen equal unit rate constants,

The positive total $T=x+y+z>0$ is conserved, so every trajectory stays on one simplex. This restriction to a positive stoichiometric class is fundamental in mass-action dynamics.⁵ Set $z=T-x-y$ to obtain the two-dimensional system

The equilibria in the triangle are $(x,y)=(0,0)$, $(T,0)$, $(0,T)$, and $(T/3,T/3)$. The first is a source, the two winner vertices are stable nodes, and the interior point is a saddle with eigenvalues $-T$ and $T/3$. One identity explains the choice:

The sign of $x-y$ is preserved and its magnitude is amplified whenever $z>0$. Positive equal initial competitors stay on the invariant diagonal and approach the saddle. The all-$Z$ vertex remains stationary. Any imbalance selects a winner. This is bistability built from elementary reactions, a conservation law, and positive invariance, without placing a Hill function into the model by hand. Related monotonicity tools extend order-based conclusions to much larger biochemical systems.⁶

Use stoichiometry, positivity and a conserved total to derive a reaction network's reduced decision geometry.

20Reuse the complete workflow

A defensible phase portrait is a reproducible argument whose biological, geometric, local, and global layers agree.

1. Function and scope. State the behavior: local recovery, two basins, a conserved family, finite-time escape, or an attracting cycle, and name the region over which it is claimed.
2. Family. State which mechanisms, signs, kinetic shapes and parameter ranges are fixed, and which are allowed to vary.
3. State. List variables, units, domain, external inputs, conserved totals, and any hidden variables removed by approximation.
4. Mechanism. Connect every addition, removal and regulatory sign to a reaction or declared effective law.
5. Global geometry. Find every scalar root or planar nullcline, label component signs, and mark invariant boundaries before integrating.
6. Local system. At every fixed point write $A=J(x^*)$ and the actual approximation $\dot z=Az$.
7. Strict test. Apply Hurwitz only when every eigenvalue has negative real part. Record positive or zero real parts separately.
8. Global certificate. Search for energy, a Lyapunov function, a conserved quantity, a trapping region, monotonicity, or stoichiometric structure.
9. Chemical audit. Check positivity, parameter independence, mass-action sign constraints and the assumptions behind every effective law.
10. Numerical audit. Compare several initial conditions, vary the time step and horizon, and preserve the code and parameters that generated every curve.⁹
11. Scope sentence. End with what has been proved, where it holds, and what remains an inference or hypothesis.

Four transfer exercises, with answers

1. A linear biochemical pair. For $A=\left[\begin{smallmatrix}-2&-3\\1&-1\end{smallmatrix}\right]$, classify the fixed point. Answer: $\tau=-3$, $\Delta=5$, and $D=-11$, so it is a stable spiral with eigenvalues $(-3\pm i\sqrt{11})/2$.

2. A semidefinite derivative. Why does $\dot V=-x^2$ not alone prove convergence to the origin in a two-state system? Answer: the derivative vanishes on the whole line $x=0$. One must determine the largest invariant subset of that line. It may be the origin, another trajectory, or a continuum of equilibria.

3. Center or cycle. Two initial conditions produce distinct closed curves that remain separated forever. Which object is supported? Answer: a center or another conserved family, not an attracting limit cycle. Perturbations do not converge back to one isolated orbit.

4. Toggle perturbation. A pulse moves a state close to the diagonal but not across it. Does the cell switch in the deterministic symmetric model? Answer: no. The diagonal is the invariant separatrix. The state returns to the attractor in its original basin. With noise, boundary crossing becomes probabilistic and requires a stochastic model.

Build a reproducible dynamical argument whose mechanism, state space, local tests, global claims and numerical checks agree.

†References

The important full texts have been inspected and archived in the project literature directory; its index records stable identifiers, retrieval dates, relevance, and checksums.