Specifying a Delay Differential Equation

From a Delayed Mechanism to an Equation

A delay differential equation (DDE) is appropriate when the rate of change of a biological state at the current time depends not only on current states, but also on one or more past state values. In a QSP model, the lag might represent a fixed period between a signal and its effect, such as the time between a precursor response and the appearance of a product. Specifying the equation well requires more than inserting \(t-\tau\): the state being sampled, the mechanism represented by the lag, the rate units, and the role of current states must all be explicit.

The goal here is to translate a stated mechanism into an equation, not to choose a particular numerical solver. Assume a deterministic model with continuous-time states and a specified finite positive delay. A fixed delay is an idealization: it says that the relevant past state is sampled at a particular lag, rather than representing a distribution of waiting times. State the mechanism and its assumptions before writing the mathematical form.

Definition: A DDE has the form \(\dot{x}(t)=f(t,x(t),x(t-\tau_1),\ldots,x(t-\tau_m),\theta,u(t))\), where \(x(t)\) is the vector of current states, \(x(t-\tau_j)\) are delayed state values, \(\tau_j>0\) are delays with units of time, \(\theta\) denotes parameters, and \(u(t)\) denotes any external inputs. Each component of \(f\) is a rate of change of the corresponding state. A valid specification defines every state, delay, and rate term, and makes the units on both sides of each equation agree.

The notation \(x(t-\tau)\) means the value of state \(x\) at an earlier time. It does not mean that the state has been shifted forward in the model, nor does it automatically mean that the state’s derivative is delayed. If a mechanism specifically depends on a past rate, that dependence must be stated and represented separately. The time argument identifies when the state is evaluated; the biological rationale identifies why that past value affects the current rate.

A specification should also make clear which states are delayed in each equation. A model may use a delayed value of one state while other terms depend on states at the current time. The state vector, parameters, and external inputs retain their usual roles: a delayed state is still a model state, and the delay is a parameter with time units. For a positive delay, specifying \(x(0)\) alone is not enough to evaluate \(x(t-\tau)\) near the start of the simulation; the required past values are called a history function and are treated in the next topic.

A Practical Formulation Sequence

Before writing symbols, describe the mechanism in words. Identify the quantity whose rate is being modeled, the event that occurs after a lag, and whether the past state drives production, inhibition, loss, or another rate. Then write a balance equation: the current rate of change is the sum of current inflows minus current outflows. A delayed mechanism changes a rate term’s argument; it does not change the fact that the derivative is evaluated at the current time.

Define the state and its units.
For each state, say what biological quantity it represents and how it is measured, such as concentration in mg/L or cell density in cells/mL.

Name the delayed mechanism.
Specify what event is delayed, the lag \(\tau\), and which past state value represents the biological signal or cohort relevant now.

Write each rate term.
Separate inflows, outflows, and feedback terms. Include current states wherever the mechanism depends on their present values.

Check units and assumptions.
Confirm that every term has units of state per time, and check that a fixed lag and the chosen functional form match the stated mechanism.

A useful dimensional rule is that \(\dot{x}(t)\) has units of state divided by time. A production term must therefore have those same units, while a first-order loss term \(kx(t)\) requires \(k\) to have units of inverse time. A delay \(\tau\) has units of time, and a ratio such as \(x(t-\tau)/K\) must be dimensionless if it appears inside a saturation function.

Worked Example: Delayed Saturable Production with Current Clearance

Consider a hypothetical mediator concentration \(C(t)\), measured in mg/L. A biological signal measured by the mediator concentration one day earlier drives saturable production, while the mediator is cleared in proportion to its current concentration. Suppose the lag is fixed, the delayed signal is nonnegative, and the clearance process is first order. Let the production capacity be \(V_{\max}=12\) mg/L/day, the half-saturation concentration be \(K=4\) mg/L, the Hill exponent be \(h=2\), and the clearance constant be \(k=0.30\) day\(^{-1}\).

The mechanism gives the equation

\frac{dC(t)}{dt} = \frac{V_{\max}}{1+\left(C(t-\tau)/K\right)^h} -kC(t), \qquad \tau=1\ \text{day}.

The numerator and denominator of the production term give units of mg/L/day, since the denominator is dimensionless. The clearance term also has units of mg/L/day because \(k\) has units day\(^{-1}\) and \(C(t)\) has units mg/L. The delayed concentration controls production; the current concentration controls clearance. This distinction is part of the biological specification, not merely a notational choice.

Suppose at a particular time \(t\), the hypothetical delayed concentration is \(C(t-\tau)=5\) mg/L and the current concentration is \(C(t)=5\) mg/L. Substitution gives

\text{production} = \frac{12}{1+(5/4)^2} = \frac{12}{1+25/16} = \frac{192}{41} \approx 4.683\ \text{mg/L/day}.

Current clearance is \(0.30(5)=1.5\) mg/L/day. Therefore the net rate is \(4.682927\ldots-1.5=3.182927\ldots\), or approximately \(3.183\) mg/L/day. The positive derivative says that concentration is increasing at this instant under the specified model and values. It does not say that concentration will continue to increase at the same rate.

Delayed Inputs to Production or Inhibition

The delayed state need not be the same quantity whose derivative is being modeled. A drug concentration, for example, can be a delayed driver of a response state. In this case the model must identify the delayed driver and the response separately. Do not insert a delay merely because the response appears slow: the assumed lag should correspond to a specific mechanism, and the form of the rate term should match that mechanism.

Worked Example: Delayed Drug-Dependent Suppression of Synthesis

Suppose a hypothetical response concentration \(R(t)\), in mg/L, is synthesized at a rate suppressed by drug concentration \(C(t)\), also in mg/L. Assume the suppression follows a simple inhibitory Emax form with maximum inhibition equal to one, and that the drug concentration affecting synthesis is sampled after a fixed delay of two days. Response is removed by first-order loss based on its current concentration. Let \(k_{\mathrm{in}}=8\) mg/L/day, \(IC_{50}=2\) mg/L, and \(k_{\mathrm{out}}=0.15\) day\(^{-1}\).

The corresponding equation is

\frac{dR(t)}{dt} = \frac{k_{\mathrm{in}}}{1+C(t-\tau)/IC_{50}} -k_{\mathrm{out}}R(t), \qquad \tau=2\ \text{days}.

The ratio \(C(t-\tau)/IC_{50}\) is dimensionless, so the synthesis term has units of mg/L/day. The loss term has the same units. Suppose \(C(t-\tau)=6\) mg/L and \(R(t)=10\) mg/L. The synthesis rate is \(8/(1+6/2)=8/4=2\) mg/L/day, and the loss rate is \(0.15(10)=1.5\) mg/L/day. Thus

\frac{dR(t)}{dt}=2-1.5=0.500\ \text{mg/L/day}.

The model predicts a positive instantaneous response rate at these specified values, even though the delayed drug level suppresses synthesis. The current response level determines the loss term; the past drug level determines the synthesis term. The equation assumes a fixed lag and this particular inhibition function. If the response were instead affected by current drug concentration, the argument in the synthesis term would be \(C(t)\), and that would be a different model.

Coupled States and Equation-by-Equation Specification

For a coupled system, define an order for the state vector and write each component equation with its own biological interpretation. The state vector order is bookkeeping; it does not imply that every component has the same delays. An equation can depend on current values of several states and on a delayed value of only one of them. This explicit mapping helps prevent indexing errors when translating a biological model into a right-hand-side function.

Worked Example: Delayed Target-Dependent Effector Recruitment

Consider a hypothetical target burden \(T(t)\) and effector-cell density \(E(t)\), both measured in cells/mL. Target cells grow logistically and are killed by current effectors. Effector recruitment responds to the target burden two days earlier, while effector cells are lost in proportion to their current density. Assume a fixed recruitment lag, nonnegative states, and the stated logistic and saturating response forms. Let \(r=0.40\) day\(^{-1}\), carrying capacity \(K=1000\) cells/mL, killing coefficient \(k=0.002\) mL/(cell·day), basal recruitment \(s=1\) cells/mL/day, maximal recruited rate \(a=6\) cells/mL/day, recruitment half-saturation \(K_T=500\) cells/mL, and effector loss \(d=0.10\) day\(^{-1}\).

Choose the state vector \(x(t)=(T(t),E(t))\). The mechanism translates to

\begin{aligned} \frac{dT(t)}{dt} &= rT(t)\left(1-\frac{T(t)}{K}\right) -kE(t)T(t),\\ \frac{dE(t)}{dt} &= s+a\frac{T(t-\tau)}{K_T+T(t-\tau)} -dE(t), \qquad \tau=2\ \text{days}. \end{aligned}

In the target equation, the logistic term and killing term each have units of cells/mL/day. In the effector equation, the saturating ratio is dimensionless, so the recruitment term has units of cells/mL/day; the loss term has those same units. The delay appears only in recruitment because the stated mechanism assigns the lag there.

For an arithmetic check, suppose \(T(t)=800\) cells/mL, \(E(t)=20\) cells/mL, and \(T(t-\tau)=500\) cells/mL. The target growth term is \(0.40(800)(1-800/1000)=64\) cells/mL/day. The killing term is \(0.002(20)(800)=32\) cells/mL/day, so \(dT/dt=64-32=32\) cells/mL/day. Effector recruitment is \(1+6(500/(500+500))=4\) cells/mL/day; effector loss is \(0.10(20)=2\) cells/mL/day. Hence \(dE/dt=4-2=2\) cells/mL/day. Both rates are positive at this instant. No inference about long-term behaviour follows from these two instantaneous rates alone.

Justifying the Placement of the Delay

A simple balance argument explains why a delayed driver appears in the rate term rather than replacing the current-time derivative. Let \(x(t)\) be a concentration produced at rate \(q\) and cleared at first-order rate \(kx(t)\). If the mechanism states that production at time \(t\) is determined by a signal \(z\) measured \(\tau\) earlier, then the amount added over a short interval \(\Delta t\) is approximately \(q(z(t-\tau))\Delta t\), while the amount removed is approximately \(kx(t)\Delta t\). The net change divided by \(\Delta t\), in the limit as \(\Delta t\) becomes small, is therefore the current derivative:

\frac{x(t+\Delta t)-x(t)}{\Delta t} \approx q(z(t-\tau))-kx(t) \quad\Longrightarrow\quad \dot{x}(t)=q(z(t-\tau))-kx(t).

This argument also shows why \(x(t-\tau)\) should not be substituted for \(x(t)\) in every term. Only the rate contribution whose mechanism depends on the past signal receives the delayed argument. Nor should \(q(z(t-\tau))\) be replaced with \(q(\dot{z}(t-\tau))\) unless the biological mechanism explicitly depends on the past rate of change of \(z\).

Assumptions, Diagnostics, and Alternatives

A fixed-delay formulation is most defensible when the relevant latency is reasonably represented by a single lag and the model’s state variables contain the information needed to define the delayed mechanism. Check that the biological description supports a fixed delay, that each delayed argument refers to the intended state, and that the rate law remains meaningful over the domain of interest. For example, a saturating function with a nonnegative concentration is naturally interpreted differently from the same expression evaluated at an impossible negative concentration.

Check dimensions term by term. A derivative and every inflow or outflow in its equation must have units of state per time. A dimensionally inconsistent equation often signals a missing rate constant or an incorrectly scaled state.
Check the time argument. Confirm that the delayed state is evaluated at \(t-\tau\), with \(\tau\) in the same time units used for \(t\). A delay entered in hours when the model time is days changes the mechanism by a factor of 24.
Check the balance structure. Production and loss terms should have signs and dependencies consistent with the biological statement. A delayed signal driving production does not automatically delay clearance.
Check limiting cases. Setting a delay to zero should recover the corresponding instantaneous-rate form. Removing a delayed mechanism (for example, setting its maximal rate to zero) should remove that contribution without changing unrelated terms.
Check the state mapping. In a vector implementation, verify that each delayed component corresponds to the named biological state rather than another state at the same index or a current-time value by mistake.

A common mistake is to write \(dC(t)/dt=P(C(t-\tau))-kC(t-\tau)\) whenever any part of the mechanism is delayed. That equation delays both production and clearance. If only the production signal is delayed and clearance acts on the present concentration, the correct structure is \(dC(t)/dt=P(C(t-\tau))-kC(t)\). The two formulations represent different mechanisms and can produce different trajectories.

An ordinary differential equation (ODE) uses current states only, so it is the simpler choice when no biologically meaningful lag is needed. A fixed-delay DDE directly represents a single specified lag. Models with added intermediate states or distributed delays can represent a spread of transit times, but they are different formulations with additional assumptions. The choice should follow the mechanism: use a fixed-delay DDE when the lag itself is the intended idealization, rather than introducing it simply to make a curve fit.

A DDE specification is not a complete simulation setup by itself. Because delayed state values are needed, the model also requires past-state information over the interval needed by its delays, as well as parameter values and any inputs. In addition, the predicted dynamics depend on the assumed fixed lag and rate laws. Sparse observations may not distinguish a delayed mechanism from other sources of sluggish response; a plausible fit alone does not establish that the delay is biologically real.

Key takeaway: In a well-specified DDE, the derivative is evaluated now, delayed arguments identify the past states that drive particular mechanisms, and every rate term has the units of its state per time. Write the biological reason for each lag explicitly and keep current-state terms current unless the mechanism says otherwise.

Check Your Understanding

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

A concentration \(A(t)\) has delayed production \(V/(1+A(t-\tau)/K)\) and current first-order clearance \(kA(t)\). Which terms use current and delayed values, and what are the units of \(V\), \(K\), and \(k\) if \(A\) is measured in mg/L?
A delayed inhibition term uses \(D(t-\tau)/IC_{50}\). What units must \(D\) and \(IC_{50}\) have for the ratio to be dimensionless?
Suppose production is delayed, but a proposed equation uses \(kA(t-\tau)\) for clearance. What mechanism does this equation imply, and how should it be corrected if clearance acts on current concentration?
For \(dR/dt=6/[1+C(t-1)/2]-0.2R(t)\), let \(C(t-1)=4\) mg/L and \(R(t)=5\) mg/L. If \(C\) and \(R\) are in mg/L and time is in days, calculate the net rate and verify the units.

Solutions. (1) Production uses \(A(t-\tau)\), while clearance uses \(A(t)\). \(V\) has units mg/L/time, \(K\) has units mg/L, and \(k\) has units inverse time. (2) Both \(D\) and \(IC_{50}\) must use the same concentration units, such as mg/L. (3) The proposed equation delays both the production driver and clearance. If clearance is current, write \(dA/dt=V/(1+A(t-\tau)/K)-kA(t)\). (4) Production is \(6/(1+4/2)=2\) mg/L/day; clearance is \(0.2(5)=1\) mg/L/day; the net rate is \(1.00\) mg/L/day. Each term has concentration-per-day units.