From simulations to possible behaviors

Build a map from biochemical structure to dynamical function: dimension decides what a vector field can do, and chemistry decides which vector fields a cell can have.

CCBS · EST 5053 · Fall 2026 · lecture 4 · core page (full) · 2026-09-10 (Thu) 10:40–12:15 · students present, lecturer-backed

You can already do something powerful. Write reactions, assemble $\dot x=\Gamma v(x)$, integrate the equations, move a slider and watch the behavior change. If that were enough, this course could end here. It is not enough, and the reason is not that simulation is bad. A simulation answers one very precise question about one declared model, parameter point and initial state.

0–5 min

Opening. Reopen a reaction-system simulation and collect what it can and cannot answer.

5–13 min

Simulation limit. Separate missing mechanisms from combinatorial explosion.

13–18 min

Goal. Draw the map from structure to function.

18–27 min

One dimension. Coordinate addition and removal, $f(x)$, a phase line and $x(t)$.

27–36 min

Memory. Construct two stable states from three crossings.

36–44 min

Bifurcation. Make tangency, hysteresis and critical slowing visible.

44–56 min

Two dimensions. Build the biological toggle before drawing its nullclines.

56–62 min

Dimension boundary. See why the plane has a short catalogue and three dimensions do not.

62–73 min

Local retreat. Write the linearized system, use Hurwitz and read the Jacobian as wiring.

73–81 min

Global retreat. Let energy settle a case where the Jacobian is silent.

81–87 min

Chemistry. Apply positivity, parameter independence and constrained kinetics.

87–90 min

Return. Fill in the structure-to-function map and read the price list forward.

90–95 min

Lecturer addendum. Connect the team's chosen network to the rest of the course.

1You can already simulate. Why is that not enough?

One trajectory is exact for the model that generated it. It cannot report a mechanism the model omitted, and it cannot exhaust a large family of parameters and initial states.

In the opening model, $x$ is protein concentration, $\mu$ is its constant production rate, $\lambda$ is its first-order removal rate, and $x_0=x(0)$ is the initial concentration. Production balances dilution or degradation at $x^{*}=\mu/\lambda$, and direct integration gives

The slider therefore changes two things together: the crossing moves as $1/\lambda$, and the recovery time changes as $1/\lambda$. The canvas is useful precisely because the mechanism and solution are already declared.

Limit 1: the discarded mechanism is not in the run

A cellular model is built by compression. Several elementary binding states become one effective Hill curve. Transcription and translation become one composite production arrow. A conserved pool may be removed from the coordinates. Each step says that some hidden process is fast, fixed or irrelevant to the chosen observable. A run cannot tell us what the discarded process would have changed, because that process is absent from the equations being run. We need claims that survive a class of omitted details.

Limit 2: the family is too large to enumerate

Ten uncertain parameters sampled at only ten settings already require $10\times10\times\cdots\times10=10^{10}$ parameter combinations, before sampling any initial condition. Finer grids multiply the problem without proving what happens between grid points. We need statements such as “every field with this sign pattern is stable” or “no scalar autonomous system can oscillate.” Those cover infinitely many numerical runs at once.

2What we actually want: structure to function

The valuable arrows connect coarse structural facts to qualitative functions. Exact parameters remain essential when the question is an exact time trace.

The left hierarchy moves from coarse structure to detail: which species interact, the signs of those interactions and conserved pools, the shapes of effective response curves, then exact parameters and the initial state. The right hierarchy moves from coarse function to detail: what state space is admissible and whether motion stays bounded, whether a steady state is stable, whether there are multiple states or a cycle, then the exact period, amplitude and trajectory.

The first promised upper arrow is concrete: whether a one-variable gene circuit can remember does not depend on knowing every rate constant. It depends first on whether a regulated addition curve, here protein production, can cross its removal curve, here dilution, three times. Sections 3 and 4 will construct that statement. Exact rate constants then decide where the crossings lie and how fast each state is approached.

3One dimension: a fixed point is a crossing

For a cellular concentration, write the net rate as a nonnegative addition term minus a nonnegative removal term. Stability becomes a comparison of slopes at their crossing.

Let $x\ge0$ be one concentration. Write its total addition rate as $f^{+}(x)\ge0$, its total removal rate as $f^{-}(x)\ge0$, and its net rate as $f(x)$:

A fixed point $x^{*}$ is an addition-removal crossing, $f^{+}(x^{*})=f^{-}(x^{*})$. Just to its right, removal must win for motion to point back; just to its left, addition must win. At a simple crossing this is the slope condition $f'(x^{*})=(f^{+})'(x^{*})-(f^{-})'(x^{*})<0$.

This split uses the physical domain of a concentration: $x\ge0$, and both accounting terms are nonnegative. A mechanical coordinate such as spring displacement is signed, and its restoring term $-kq$ need not be separated into nonnegative addition and removal. The distinction is structural, not cosmetic. Section 10 returns to the restrictions imposed by chemical state variables.

1. Plot $f(x)$, or plot $f^{+}(x)$ and $f^{-}(x)$ and compare them.
2. Mark every zero of $f$.
3. On each interval between zeros, draw motion right for $f>0$ and left for $f<0$.
4. Translate the phase line into time traces from several initial conditions. The traces must approach, depart or escape exactly as the arrows predict.

For $\dot x=-x$, the field crosses with negative slope, both arrows point to zero and every trajectory decays exponentially. For $\dot x=x^2$, the field only touches zero. The origin attracts from the left and repels from the right, so $f'(0)=0$ cannot classify it. Separating variables gives

When $x_0>0$, escape occurs at $t=1/x_0$. The phase line got the direction right but carried no clock. This signed mathematical example is not yet a closed concentration model; section 10 will make that boundary explicit.

The third row has two attracting roots separated by one repelling threshold. It also proves the first impossibility result. Between consecutive roots the sign of $f$ cannot change, so a nonstationary scalar trajectory is monotone there. To return to an earlier value it would have to reverse direction, which requires crossing a zero, but a unique trajectory cannot cross a fixed point. A scalar autonomous system can switch, but it cannot sustain a nonconstant cycle.¹

This is the first structural arrow. Deform the addition or removal curve without changing the number, order and transverse character of their crossings, and the same phase-line portrait survives. The exact trajectories change; the qualitative basins do not.

4To remember, cross three times

Memory is a construction problem. Two persistent states require stable, unstable, stable crossings, so the regulated addition curve must locally outrun removal and later saturate.

Consider a positively autoregulated gene whose protein activates its own production. Fast redistribution among promoter states is compressed into an activating Hill fraction, while protein removal is dominated by first-order dilution. In dimensionless concentration $x$ and dimensionless time $\tau$, one concrete teaching model is

The Hill term is a lumped effective law. It assumes promoter binding or redistribution is fast relative to protein accumulation, and the exponent represents the chosen cooperative occupancy model. It is not an elementary four-molecule collision and it is not a universal law for activation.

The design logic does not begin with those four numbers. It begins with shape. Basal production keeps the low state away from zero. Positive feedback makes production rise. Sufficient effective cooperativity lets its slope exceed the removal slope over an interval. Saturation makes it flatten so removal can win again. In this construction, memory is bought with a sufficiently steep cooperative response. Other molecular mechanisms can create the required sigmoid, so cooperativity is a route to the shape rather than a proof of one unique mechanism.

5When two curves touch, a switch is born

A saddle-node bifurcation is a tangency in the addition-removal picture. Draw the vector field below, at and above tangency before stacking those cases into branches.

Let $r$ be a control parameter and let $x$ be a signed local displacement from the collision point. The saddle-node normal form is

In the underlying biological picture, two curves touch when both their values and slopes agree:

For $r>0$, the roots are $x^{*}=\pm\sqrt r$. Since $f_x=-2x$, the positive root is stable and the negative root is unstable. At $r=0$ they collide, and for $r<0$ neither remains. The fold creates one stable-threshold pair; when another stable branch already exists elsewhere, that event opens or closes a bistable switching range. A bifurcation is this change in the organization of all trajectories, not merely a bend in one plotted path.

Two experimental signatures

Hysteresis. A full sigmoidal switch commonly has two folds. Increasing an input destroys one state at one fold; decreasing it destroys the other state at the other fold. The switching threshold therefore depends on direction and history.

Critical slowing. Near the stable root, let $u=x-\sqrt r$ be a small displacement. Its linear term is $\dot u=-2\sqrt r\,u$, hence

Moving from $r=1$ to $r=0.01$ changes the relaxation time from 0.5 to 5. A hundredfold approach to the fold gives a tenfold slowing in this normal form. That timing prediction can be tested, provided noise, sampling and the valid parameter range are checked.²

6A second variable buys going around

In a plane, trajectories can turn around fixed points and a saddle can separate basins. Nullclines make those motions readable as addition-removal balances for individual species.

The nullcline of species $i$ is the set where its total addition equals its total removal, so its own derivative is zero there. A trajectory crosses that curve using the other component. Only an intersection of all nullclines is a fixed point. This is the right place to build a biological two-variable switch.

Let gene $G_X$ produce repressor $X$ and gene $G_Y$ produce repressor $Y$. Protein $X$ represses production from $G_Y$; protein $Y$ represses production from $G_X$. Both proteins are lost through molecular degradation, growth dilution or both. At the composite level,

For the arm producing $X$, let $G_X$ be free promoter, $C_{G_XY}$ its complex with repressor $Y$, $q_{G_X}$ the conserved promoter total, $K_Y$ an effective repression scale and $n_Y$ an effective cooperativity. A rapid-equilibrium occupancy model writes

Solving for the free promoter fraction gives

If promoter redistribution is fast relative to expression and composite production is proportional to the free promoter fraction, then the slow protein concentrations obey

Here $\alpha_X,\alpha_Y$ are maximal production rates; $K_X,K_Y$ are repression scales; $n_X,n_Y$ are effective cooperativities; and $\delta_X,\delta_Y$ combine degradation and dilution. Equation (11) is not elementary mass action. It is a reduced law for a specific biological scenario, and its assumptions stay attached.

Scale concentrations by repression constants and time by a common removal time. A symmetric teaching case is

1. Set $\dot x=0$: $x=4/(1+y^2)$. Horizontal motion pauses on this $x$-nullcline.
2. Set $\dot y=0$: $y=4/(1+x^2)$. Vertical motion pauses on this $y$-nullcline.
3. At every intersection, both components vanish. Testing one point in each region supplies the arrow directions.
4. The two outer intersections attract. The middle intersection is a saddle, and its stable manifold is the basin boundary. Symmetry makes the diagonal $x=y$ invariant here.

A mutual-repression toggle was built in Escherichia coli using two repressors and two chemical inputs to set the state.³ Equation (12) is a teaching reduction, not a fit to that device. What transfers is the dynamical claim: two attractors, their basins and a separating threshold make history matter.⁴

7A curve packs freely in three dimensions

Trajectories cannot cross, but the force of that rule depends on dimension. A line orders motion; a closed curve separates a plane; a curve in space can pass around itself.

This packing argument explains why one-dimensional behavior is almost completely read from signs and why generic planar fixed points and bounded recurrent behavior have a short qualitative catalogue: nodes, saddles, spirals, centers and isolated cycles. In a bounded planar region, the Poincaré–Bendixson theorem can turn “no fixed point remains in the limiting set” into the existence of a periodic orbit. Its hypotheses concern a two-dimensional invariant set; they do not automatically apply to a three-species network.

A three-repressor ring illustrates the new room:

Here $p_i$ is the dimensionless level of repressor $i$, $\alpha$ is maximal production and $n$ is effective cooperativity. This cyclic wiring was realized as the repressilator, and later designs improved the persistence and synchronization of the oscillation.⁵⁶ Three dimensions also permit chaos, but permission is not a mechanism: lack of structure is not a structure. When general geometry becomes permissive, we retreat to local structure and global certificates.

8Close to a fixed point, every network is linear

The Jacobian is not merely a table of derivatives. Evaluated at one fixed point, it is the matrix of a new first-order dynamical system for small displacements.

Lecture 3 ended with the cell as a crowded chemical plant and with a question that balance alone could not answer: after a perturbation, does a steady state recover? Stability is the first requirement for a persistent operating point and the local prerequisite for discussing thresholds, switches and cycles. The first retreat from an unmanageable model family is therefore to keep only what every smooth network looks like near one fixed point.

Let the state vector $x$ obey $\dot x=F(x)$, and let $x^{*}$ be a fixed point with $F(x^{*})=0$. Define the displacement vector $u=x-x^{*}$. Taylor expansion gives

The matrix entry $A_{ij}$ is the derivative of component $F_i$ with respect to state $x_j$, evaluated at the fixed point. The canonical linearized system discards the higher-order remainder and uses a new displacement vector $z$:

Equation (15), not the derivative array alone, is the linearization. It describes nearby displacement from one fixed point. A different fixed point generally produces a different matrix.

If $Aw=\lambda w$, motion along eigenvector $w$ has scalar amplitude $a(t)$ satisfying $\dot a=\lambda a$. Thus

The real part of $\lambda_j$ sets the exponential envelope; its imaginary part sets rotation. Negative real part shrinks a mode, positive real part grows it. Every small displacement is assembled from the modes, so all of them must shrink for the fixed point to recover from every small perturbation.

For a real two-by-two matrix, let $\tau=\operatorname{tr}A$, $\Delta=\det A$ and $D=\tau^2-4\Delta$. Since the eigenvalues sum to $\tau$ and multiply to $\Delta$,

The intuition is direct. A negative determinant means two real eigenvalues of opposite sign, so one direction grows. A positive determinant puts both real eigenvalues on the same side, or gives a conjugate pair. A negative trace then places their sum, and therefore both real parts, on the decaying side.

For the toggle in equation (12),

Growth dilution contributes a term $-\lambda_i x_i$ to a concentration equation and therefore contributes $-\lambda_i$ to the corresponding diagonal entry before other interactions are added. It supplies local self-removal, but it is not a universal stability guarantee. In the toggle, the negative diagonal is protein removal and the two negative off-diagonal entries are mutual repression. In any biochemical model, $A_{ij}$ asks how species $j$ changes species $i$'s net rate. The sign pattern of the Jacobian is the local wiring diagram.

9When the linear test says nothing

Identical purely imaginary Jacobians can hide attraction, neutral cycles or repulsion. A monotone scalar function can settle what the discarded nonlinear terms do.

Compare three systems written in polar radius $r$ and angle $\theta$:

At the origin, all three have Jacobian $\left[\begin{smallmatrix}0&-1\\1&0\end{smallmatrix}\right]$ and eigenvalues $\pm i$. The cubic radial term is invisible at first order and decides three different fates.

The spring now enters with two jobs in the same narrative: it explains what a state must contain, and its energy supplies a global certificate. Let $q$ be displacement, $m$ mass, $k$ spring stiffness and $c$ viscous drag. Force balance gives

Position alone is not a state because the mass may pass the same $q$ moving left or right. Set $v=\dot q$ and let the state vector be $z=(q,v)$. The second-order equation becomes an autonomous first-order system:

This conversion is exact, not a linearization. More generally, a higher-order law or a process with memory becomes first order by enlarging the state until the present determines the instantaneous future. Hidden promoter states and delayed intermediates pose the same modeling question.

Define the mechanical energy

Differentiate along the full dynamics:

With $c=0$, energy is conserved and each nonzero initial condition stays on its own ellipse. The origin is stable but not attracting. With $c>0$, energy never increases. The derivative is zero on the entire line $v=0$, so strict decrease cannot be claimed. But at any point on that line with $q\ne0$, $\dot v=-(k/m)q\ne0$; the trajectory immediately leaves. Only the origin remains invariant there, so LaSalle's invariance argument gives convergence.

10Chemistry forbids most vector fields

The spring was useful borrowed physics, but its signed rotation is not a closed concentration model. Mass-action networks inherit positivity, independent parameters and constrained polynomial terms.

Write elementary reaction $i$ as

Here $X_j$ is species $j$, $x_j$ is its concentration, $\alpha_{ji}$ and $\beta_{ji}$ are reactant and product stoichiometric coefficients, $k_i$ is the rate constant and $v_i$ is the elementary mass-action flux. If reaction $i$ consumes $X_j$, then $\alpha_{ji}\ge1$, so its flux contains $x_j$ and vanishes on the boundary $x_j=0$. Consumption therefore cannot push a concentration negative.

Three restrictions to carry forward

Positivity. Every face of the nonnegative orthant must point inward or tangent. The scalar law $\dot x=x-\mu$ with $\mu>0$ and the pure rotation in Figure 10 both fail at the boundary.

Independent parameters. Distinct elementary reactions have distinct rate constants. A behavior that requires two unrelated constants to match exactly is tuned, not structural. A dual-rail representation can encode a signed variable as $q=x_+-x_-$, but an exact arbitrary linear system generally requires matched rates in the two rails. Biology may regulate those rates, but the equality then needs a mechanism of its own.

Constrained kinetic form. Elementary mass action produces polynomials with signs tied to what each reaction consumes and produces. Effective non-polynomial laws such as Hill and Michaelis-Menten functions arise by eliminating fast binding or catalytic states. Timescale separation is therefore a genuine route past the polynomial wall, not decorative curve fitting. Lecture 6 will ask exactly when that elimination preserves the observable we care about.

These restrictions support high-dimensional conclusions that a generic vector field cannot offer. Complex-balanced mass-action networks admit strong Lyapunov structure, while signed interaction graphs can reveal monotone or near-monotone organization.⁹¹⁰

11The price list

Analysis did not replace simulation. It identified which questions can be settled before exact parameters are known and which still require the bottom-level model.

The cell-as-chemical-plant question now has a more precise answer. Positivity keeps every concentration in its admissible region. Under a constant-growth approximation, dilution contributes self-limiting diagonal terms. Saturation caps many effective production fluxes. Conserved pools confine motion to lower-dimensional slices. None of these alone is a universal stability theorem, and cells deliberately switch and oscillate. Together they explain why unbounded or stalled behavior is not the default, and why interesting behavior must be built.

A student team can start here

1. State one biological function: recovery, memory, commitment or oscillation.
2. Draw the reaction or mechanism network and name every approximation used to compress it.
3. Declare the state, physical domain, conserved totals, parameters and initial conditions.
4. Write the vector field and connect every addition, removal and regulatory sign to the network.
5. Build the global picture first: scalar crossings or planar nullclines, fixed points, signs and candidate basins.
6. At each fixed point, write $A=J(x^{*})$ and the system $\dot z=Az$; then state exactly what Hurwitz does or does not prove.
7. Search for positivity, conservation, trapping or a Lyapunov function when the claim is global or the Jacobian is silent.
8. Use simulation last to test timing, basin boundaries, parameter variation and the assumptions that produced the effective law.

The forward price list is now visible. Memory and commitment lead to multistability and bifurcations. A clock leads to genetic and metabolic oscillators. Return after a perturbation leads to adaptation. Wiring-level conclusions lead to reaction-order geometry. Timescale separation explains when the effective Hill laws used here can be trusted.

Simulation tells you what one cell did once. Structure tells you what every cell of that design can and cannot do, and that is the only kind of answer an experiment can be designed against.

†References

These sources support the biological examples and the mathematical claims used in the lecture. Important full texts have been inspected and archived with the project.