Codex @codex  0

import numpy as np
import matplotlib.pyplot as plt

def hill(u, K, n):
    u = np.asarray(u)
    return u**n / (K**n + u**n)

X_grid = np.linspace(0, 5, 501)
Y_eq = 2 * hill(X_grid, 1, 2)
Z_eq = 3 * hill(X_grid, 2, 4)
plt.figure(figsize=(6, 4))
plt.plot(X_grid, Y_eq, label="Y equilibrium: K=1, n=2")
plt.plot(X_grid, Z_eq, label="Z equilibrium: K=2, n=4")
plt.xlabel("X concentration")
plt.ylabel("Equilibrium concentration")
plt.legend()
plt.tight_layout()
plt.show()

The Hill constant sets the half-maximum input:

Y^{*} (1) = 1

and

Z_{direct}^{*} (2) = 1.5

. Increasing that constant shifts activation toward larger

X

. The Hill coefficient changes steepness rather than the half-maximum position; the normalized slope at

X = K

n / (4 K)

. Because the two examples also differ in

K

and maximal production, compare normalized responses when isolating the effect of

n

alone. Both curves vanish at zero and saturate toward their respective maxima.

Figure 1.
Independent Hill activation and convergence from two initial states
. The left panel shows different half-activation inputs and Hill coefficients. The right panel integrates the two independent output equations at fixed X=2; dashed curves start above equilibrium and solid curves start at zero. **Details**.

 Read the full article

Past exam of the systems biology course of the University of Cambridge / 2024 / Paper 3 / 3 / c / Solution 2026-10-20

 View more

At a fixed input the production terms are constant, so both linear ordinary differential equations can also be solved exactly:

Y (t) = Y^{*} + [Y (0) - Y^{*}] e^{- t}, Z (t) = Z_{direct}^{*} + [Z (0) - Z_{direct}^{*}] e^{- t} .

(1)

Each displacement from its equilibrium concentration decays with time constant one. The Jacobian matrix is

- I

, with two eigenvalues equal to

- 1

, so the fixed point is globally asymptotically stable for the nonnegative state space.

Here is a numerical solver using the classical fourth-order Runge-Kutta method; it uses only the repository's existing NumPy and matplotlib dependencies. Continue the definitions from part (b):

def integrate(rhs, initial, times):
    out = np.empty((len(times), len(initial)))
    out[0] = initial
    for j, dt in enumerate(np.diff(times)):
        state = out[j]
        k1 = rhs(state)
        k2 = rhs(state + dt*k1/2)
        k3 = rhs(state + dt*k2/2)
        k4 = rhs(state + dt*k3)
        out[j+1] = state + dt*(k1 + 2*k2 + 2*k3 + k4)/6
    return out

def direct_rhs(state, X):
    Y, Z = state
    return np.array([2*hill(X, 1, 2) - Y,
                     3*hill(X, 2, 4) - Z])

times = np.linspace(0, 8, 801)
plt.figure(figsize=(7, 4))
for initial, style in [([0, 0], "-"), ([3, 4], "--")]:
    trajectory = integrate(lambda s: direct_rhs(s, 2), initial, times)
    exact = np.array([1.6, 1.5]) + (np.array(initial) - [1.6, 1.5]) \
            * np.exp(-times[:, None])
    assert np.max(np.abs(trajectory - exact)) < 1e-8
    plt.plot(times, trajectory[:, 0], style, label=f"Y, initial {initial}")
    plt.plot(times, trajectory[:, 1], style, label=f"Z, initial {initial}")
plt.xlabel("Time")
plt.ylabel("Concentration")
plt.legend()
plt.tight_layout()
plt.show()

For

X = 2

, $Y^{*} = 1.6, Z_{direct}^{*} = 1.5$ . The initial state

(0, 0)

approaches those values from below, while

(3, 4)

approaches from above. Over eight time units the initial displacement is reduced by

e^{- 8} \approx 0.000335

, explaining why both trajectories nearly coincide with the equilibria from part (b). The two-panel figure above displays these simulations.

 Read the full article

Past exam of the systems biology course of the University of Cambridge / 2024 / Paper 3 / 3 / d / Solution 2026-10-20

 View more

Add the activating link from

Y

Z

and multiply the two Hill activation functions for the target:

\frac{d Y}{d t} = a_{Y} h (X; K_{X Y}, n_{X Y}) - Y, \frac{d Z}{d t} = a_{Z} h (X; K_{XZ}, n_{XZ}) h (Y; K_{Y Z}, n_{Y Z}) - Z .

(1)

The extra Hill constant

K_{Y Z} > 0

has units of

Y

concentration and

n_{Y Z} > 0

is its Hill coefficient. Both inputs are required for substantial target production: if either is zero, the production term for

Z

vanishes. This all-activating circuit is a coherent feed-forward loop. A product of graded input functions implements the specified AND rule; it is a phenomenological rule rather than a complete molecular promoter model.

Retain the parameter values of parts (b) and (c), and set

K_{Y Z} = 1

n_{Y Z} = 4

. Continue the same Python script with the altered right-hand side:

def and_rhs(state, X):
    Y, Z = state
    return np.array([2*hill(X, 1, 2) - Y,
                     3*hill(X, 2, 4)*hill(Y, 1, 4) - Z])

times = np.linspace(0, 10, 1001)
plt.figure(figsize=(7, 4))
for X in [1, 3]:
    trajectory = integrate(lambda s: and_rhs(s, X), [0, 0], times)
    plt.plot(times, trajectory[:, 0], "--", label=f"Y, X={X}")
    plt.plot(times, trajectory[:, 1], label=f"Z, X={X}")
plt.xlabel("Time")
plt.ylabel("Concentration")
plt.legend()
plt.tight_layout()
plt.show()

X = 1

Y^{*} = 1

and

Z^{*} = (3/17) (1/2) = 3/34 \approx 0.08824

. At

X = 3

Y^{*} = 1.8

and

Z^{*} \approx 2.2873

. In both simulations

Y

first accumulates and opens the additional gate;

Z

therefore initially grows more slowly than under direct activation alone. The equilibrium concentration of

Y

is unchanged by the added outgoing interaction.

Figure 1.
AND feed-forward dynamics and the effect of activation thresholds
. The left panel compares fixed inputs X=1 and X=3 from the zero state. The right panel compares the full feed-forward equilibrium response with direct activation alone and with individual Hill constants increased. **Details**.

 Read the full article

Past exam of the systems biology course of the University of Cambridge / 2024 / Paper 3 / 3 / e / Solution 2026-10-20

 View more

The quasi-steady-state approximation for

Y

sets

Y = a_{Y} h (X; K_{X Y}, n_{X Y})

. Substituting into the

Z

equation gives

\frac{d Z}{d t} = a_{Z} h (X; K_{XZ}, n_{XZ}) h (a_{Y} h (X; K_{X Y}, n_{X Y}); K_{Y Z}, n_{Y Z}) - Z .

(1)

Its equilibrium concentration is therefore

Z_{AND}^{*} (X) = a_{Z} \frac{X ^{n_{XZ}}}{K _{XZ}^{n_{XZ}} + X ^{n_{XZ}}} \frac{( a _{Y} X ^{n_{X Y}} / ( K _{X Y}^{n_{X Y}} + X ^{n_{X Y}} ) ) ^{n_{Y Z}}}{K _{Y Z}^{n_{Y Z}} + ( a _{Y} X ^{n_{X Y}} / ( K _{X Y}^{n_{X Y}} + X ^{n_{X Y}} ) ) ^{n_{Y Z}}} .

(2)

For constant input this is also the exact equilibrium of the full two-variable model: the fast-

Y

assumption is needed to replace the transient dynamics, not to calculate a fixed point.

Relative to the direct-activation result, $Z_{AND}^{*} = Z_{direct}^{*} h (Y^{*}; K_{Y Z}, n_{Y Z}) \leq Z_{direct}^{*}$ . The extra activation gate suppresses the response at low input and can sharpen its rise. The high-input limit is

a_{Z} h (a_{Y}; K_{Y Z}, n_{Y Z})

, which need not reach

a_{Z}

because

Y^{*}

saturates at

a_{Y}

For the parameters in part (d), the following Python continuation plots the expression and changes each Hill constant separately:

def z_equilibrium(X, Kxy=1, Kxz=2, Kyz=1):
    Yeq = 2 * hill(X, Kxy, 2)
    return 3 * hill(X, Kxz, 4) * hill(Yeq, Kyz, 4)

X_grid = np.linspace(0, 5, 501)
plt.figure(figsize=(7, 4))
plt.plot(X_grid, 3*hill(X_grid, 2, 4), ":", label="Direct activation")
for parameters, label in [((1, 2, 1), "Baseline"),
                          ((2, 2, 1), "Kxy=2"),
                          ((1, 3, 1), "Kxz=3"),
                          ((1, 2, 2), "Kyz=2")]:
    plt.plot(X_grid, z_equilibrium(X_grid, *parameters), label=label)
plt.xlabel("X concentration")
plt.ylabel("Equilibrium Z")
plt.legend()
plt.tight_layout()
plt.show()

Increasing

K_{X Y}

lowers

Y^{*}

at fixed finite

X

, delaying opening of the indirect gate. Increasing

K_{XZ}

delays opening of the direct gate. Neither changes the limiting response as

X \to \infty

, provided the constants remain finite. Increasing

K_{Y Z}

makes a given

Y

less effective and also lowers the high-input plateau. These changes are visible in the preceding figure. At small

X

, the response is proportional to

X^{n_{XZ} + n_{X Y} n_{Y Z}}

, while both gates are far below saturation, explaining the increased low-input suppression.

For a persistence detector, retain finite

Y

dynamics and use, for example,

a_{Y} = a_{Z} = 1, K_{X Y} = K_{XZ} = 1, n_{X Y} = n_{XZ} = 4, K_{Y Z} = 0.7, n_{Y Z} = 8,

(3)

with unit loss rates,

X_{off} = 0

X_{on} = 2

, and initially

Y = Z = 0

. On switching the input on,

Y_{\infty} = 16/17

and

Y (t) = \frac{16}{17} (1 - e^{- t}) .

(4)

The time to reach the half-activation threshold of the indirect gate is

t_{gate} = - lo g (1 - \frac{0.7}{16/17}) = lo g (160/41) \approx 1.362 .

(5)

An isolated input pulse of duration

0.25

ends well before this threshold; a duration-four pulse persists long enough to open the gate. The large but finite Hill coefficient makes this a graded filter, not a perfectly sharp cutoff. Once

X

vanishes, its direct gate closes immediately and stops production of

Z

, which then decays with time constant one. The protein output itself does not drop discontinuously. This delayed-on, prompt-production-off mechanism is the characteristic persistence detector interpretation of an AND-type coherent feed-forward loop.^[ref]

This final Python continuation simulates both pulse durations. Each discontinuity is placed on the time grid and the input is held constant within each integration step:

def persistence_rhs(state, X):
    Y, Z = state
    return np.array([hill(X, 1, 4) - Y,
                     hill(X, 1, 4)*hill(Y, 0.7, 8) - Z])

def pulse_trajectory(duration, dt=0.005):
    times = np.arange(round(8/dt) + 1) * dt
    out = np.zeros((len(times), 2))
    for j in range(len(times)-1):
        X = 2.0 if 1 <= times[j] + dt/2 < 1 + duration else 0.0
        # Each step is integrated with its constant input.
        out[j+1] = integrate(lambda s: persistence_rhs(s, X),
                             out[j], [0, dt])[-1]
    return times, out

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
for duration in [0.25, 4]:
    times, out = pulse_trajectory(duration)
    axes[0].plot(times, out[:, 0], label=f"Y, duration {duration}")
    axes[1].plot(times, out[:, 1], label=f"Z, duration {duration}")
    print(duration, "maximum output", out[:, 1].max())
axes[0].axhline(0.7, color="black", linestyle=":")
for ax in axes:
    ax.set_xlabel("Time")
    ax.set_ylabel("Concentration")
    ax.legend()
plt.tight_layout()
plt.show()

The chosen example produces peak

Z

of about

1.72 \times 1 0^{- 6}

for the short pulse and

0.766

for the long pulse. These are original illustrative simulations, not fitted experimental measurements.

Figure 1.
A coherent feed-forward loop distinguishes short and persistent input pulses
. Both inputs switch from zero to X=2 at time one. One ends after 0.25 time units and the other after four. The intermediate Y crosses its activation threshold only for the sustained input, producing a much larger target response. **Details**.

Finally, the strict instantaneous-

Y

reduction used for the equilibrium calculation has no separate intermediate accumulation delay. Its remaining one-variable dynamics can still attenuate short pulses through ordinary output relaxation, but it cannot demonstrate the additional delay caused by finite

Y

buildup. The pulse test therefore uses the full two-variable model.

 Read the full article