17  Intervening in Deterministic Systems

Status: Draft

v0.4

17.1 Introduction

This chapter shows how to intervene in deterministic dynamical systems, modifying ODEs through the \(do(\cdot)\) operator and evaluating policies. This is part of Doing in the Dynamical stratum, how can we intervene/modify dynamic processes? In continuous-time causal models (Peters et al.’s causal kinetic models), an intervention replaces one or more coordinate mechanisms in the governing equations (hard point interventions, soft mechanism edits, or structural removals), then integrates the modified system forward (Peters et al. 2022).

17.2 Interventions in ODEs

17.2.1 The Intervention Operator

When we intervene \(do(A_t = a)\) in a deterministic system, we modify the structural equation for \(A_t\):

\[ A_t \coloneqq a \quad \text{(instead of } \pi(H_t, U^A_t) \text{)} \]

This structural change propagates through the ODE system.

17.2.2 How Interventions Modify System Dynamics

Steps:

  1. Set intervention: \(do(A_t = a)\) for all \(t\) in intervention period
  2. Modify ODE: Replace \(A_t\) with constant \(a\) in the ODE
  3. Solve ODE: Integrate the modified system forward
  4. Observe effects: How does the system evolve under intervention?

17.2.3 Example: Vaccination in SIR Model

Consider an SIR model with vaccination intervention:

\[ \begin{aligned} \frac{dS}{dt} &= -\beta S I - \nu \cdot do(V = 1) \\ \frac{dI}{dt} &= \beta S I - \gamma I \\ \frac{dR}{dt} &= \gamma I + \nu \cdot do(V = 1) \end{aligned} \]

Under intervention \(do(V = 1)\), vaccination rate \(\nu\) is applied, moving individuals directly from \(S\) to \(R\).

17.2.4 Implementation: Intervention in ODEs

We can implement interventions in ODEs using DifferentialEquations.jl. Here’s an example with the SIR model:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using OrdinaryDiffEq CairoMakie

# SIR model parameters
β = 0.3  # Transmission rate
γ = 0.1  # Recovery rate
ν = 0.05  # Vaccination rate (when intervention is active)

# Original SIR model (no intervention)
function sir_original!(du, u, p, t)
    """SIR epidemic model: Susceptible → Infected → Recovered."""
    S, I, R = u
    du[1] = -β * S * I  # dS/dt
    du[2] = β * S * I - γ * I  # dI/dt
    du[3] = γ * I  # dR/dt
end

# SIR model under intervention do(V = 1) (vaccination active)
function sir_intervened!(du, u, p, t)
    """SIR model with vaccination intervention: do(V = 1)."""
    S, I, R = u
    # Intervention: do(V = 1) means vaccination is active
    du[1] = -β * S * I - ν  # dS/dt: vaccination moves S → R
    du[2] = β * S * I - γ * I  # dI/dt: unchanged
    du[3] = γ * I + ν  # dR/dt: vaccination adds to recovered
end

# Initial conditions
u0 = [0.99, 0.01, 0.0]  # 99% susceptible, 1% infected, 0% recovered
tspan = (0.0, 100.0)

# Solve original system
prob_original = ODEProblem(sir_original!, u0, tspan)
sol_original = solve(prob_original, Tsit5())

# Solve intervened system
prob_intervened = ODEProblem(sir_intervened!, u0, tspan)
sol_intervened = solve(prob_intervened, Tsit5())

# Compare peak infection
peak_original = maximum([u[2] for u in sol_original.u])
peak_intervened = maximum([u[2] for u in sol_intervened.u])
println("Peak infection rate:")
println("  Original: ", round(peak_original, digits=3))
println("  Intervened (do(V=1)): ", round(peak_intervened, digits=3))
println("  Reduction: ", round((1 - peak_intervened/peak_original) * 100, digits=1), "%")
Peak infection rate:
  Original: 0.302
  Intervened (do(V=1)): 0.033
  Reduction: 89.2%
Figure 17.1: Intervention in SIR model: comparing original and intervened trajectories

17.2.5 CausalDynamics continuous CDM (do_pin)

The vaccination example above edits the RHS by hand. CausalDynamics provides the same continuous-CDM vocabulary on a named spec: hard pin (do_pin), initial-condition (do_ic), soft force (do_force), and RHS replacement (do_rhs), with SciML-native pin callbacks (Peters et al. 2022).

include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using CausalDynamics OrdinaryDiffEq Graphs

spec = ContinuousCDMSpec(
    [:prey, :predator];
    parents = Dict(
        :prey => [:prey, :predator],
        :predator => [:prey, :predator],
    ),
)
function lotka!(du, u, p, t)
    X, Y = u
    du[1] = p.r * X - p.α * X * Y
    du[2] = p.β * X * Y - p.δ * Y
    return nothing
end
p = (r = 1.0, α = 0.1, β = 0.02, δ = 0.5)
sol = solve_cdm(spec, lotka!, [40.0, 9.0], (0.0, 8.0), p)
sol_pin = solve_cdm(
    spec, lotka!, [40.0, 9.0], (0.0, 8.0), p;
    intervention = do_pin(:predator, 5.0),
)
g_cdm = continuous_cdm_graph(spec)
(
    terminal_state(spec, sol).predator,
    terminal_state(spec, sol_pin).predator,
    ne(g_cdm),
)
(7.0788894614717615, 5.0, 2)

Continuous CDM Lotka–Volterra: observational trajectory versus hard pin do_pin(:predator, 5). Predator stays flat under the pin; prey responds through the remaining mechanism.

17.2.6 Inferring ODE parents across environments

When the parent graph is unknown but several experimental environments are observed, infer_ode_parents ranks candidate parents by leave-one-environment stability of a linear-in-parameters mechanism for \dot Y (CausalKinetiX reference method (Pfister et al. 2019)). The ranking maps onto the same ContinuousCDMSpec parent field used above.

include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using CausalDynamics Random Graphs

rng = Random.Xoshiro(1)
L = 40
times = collect(range(0.0, 4.0; length = L))
trajectories = Matrix{Float64}[]
env = Int[]
α, β = -0.5, 1.2
for e in 1:4
    for _ in 1:2
        X1 = 0.5 + 0.3 * e .+ 0.1 .* sin.(times .+ 0.2 * e) .+ 0.02 .* randn(rng, L)
        X2 = randn(rng, L)
        Y = zeros(L)
        Y[1] = 0.2 + 0.1 * e
        dt = times[2] - times[1]
        for ℓ in 1:(L - 1)
            dY = α * Y[ℓ] + β * X1[ℓ] + 0.02 * randn(rng)
            Y[ℓ + 1] = Y[ℓ] + dt * dY
        end
        push!(trajectories, vcat(Y', X1', X2'))
        push!(env, e)
    end
end

ranking = infer_ode_parents(times, trajectories, env, 1; max_size = 2, K = 3)
spec_ode = ode_parent_ranking_to_continuous_spec(ranking, [:Y, :X1, :X2]; max_parents = 1)
g_ode = continuous_cdm_graph(spec_ode)
(ranking.variable_scores, spec_ode.parents[:Y], ne(g_ode))
([0.3333333333333333, 1.0, 0.3333333333333333], [:X1], 1)

ODE parent ranking across environments. Left: multi-environment trajectories (true driver X_1 and spurious X_2). Centre: inclusion scores (higher is better). Right: recovered continuous-CDM parent graph X_1 \rightarrow Y.

17.3 Taxonomy of Continuous-Time Interventions

The interventions introduced above can be formalised into a taxonomy that distinguishes four fundamental types. Each type has distinct mathematical properties and appropriate domains of application.

17.3.1 Point Interventions

Point interventions are instantaneous interventions at a single time \(t^*\):

\[ X(t^*) = x^* \]

Equivalently, they can be represented as Dirac delta forcing in the dynamics:

\[ \dot{X} = f(X) + c \cdot \delta(t - t^*) \]

where \(\delta(t - t^*)\) is the Dirac delta distribution concentrated at \(t^*\). The constant \(c\) determines the magnitude of the instantaneous impulse.

Mathematical properties:

  • The intervention creates a discontinuity: the state jumps from its pre-intervention value to \(x^*\) (or receives an impulse)
  • The modified ODE is solved in two phases: integrate from initial conditions to \(t^*\), apply the intervention, then integrate forward from the new state
  • Point interventions apply when the intervention duration is negligible compared to the system’s characteristic timescales

Biological example: A bolus drug injection delivers a fixed dose instantaneously. The drug concentration jumps at the injection time, then decays according to the system’s natural dynamics (e.g., clearance).

17.3.2 Sustained Interventions

Sustained interventions hold a variable at a constant value over an interval \([t_1, t_2]\):

\[ X(t) = x^* \quad \text{for } t \in [t_1, t_2] \]

This is equivalent to clamping: the variable is removed from the dynamics during the intervention period. Downstream variables that depend on \(X\) receive the fixed value \(x^*\) instead of the natural evolution of \(X\).

Mathematical formulation: The modified ODE has boundary conditions \(X(t_1) = x^*\) and \(X(t_2) = x^*\), and \(X\) is treated as exogenous (not governed by \(\dot{X} = f(X)\)) during \([t_1, t_2]\). The system reduces to fewer state variables during the intervention.

Biological example: Continuous drug infusion maintains a constant drug concentration (or infusion rate) over time. The intervention sustains the treatment level rather than delivering it as a single pulse.

17.3.3 Stochastic Interventions

Stochastic interventions specify that the intervened variable is drawn from a distribution \(G\) rather than fixed at a value:

\[ do(X \sim G) \quad \text{instead of } do(X = x) \]

This generalises the standard \(do(X = x)\) to allow for uncertainty or heterogeneity in the intervention.

Connection to epidemiology: In shift interventions, the distribution of \(X\) is shifted (e.g., \(do(X \sim X + \delta)\)). Stochastic interventions formalise this by specifying a target distribution \(G\).

Mathematical formulation:

\[ P(Y \mid do(X \sim G)) = \int P(Y \mid do(X = x)) \, dG(x) \]

The interventional distribution under \(do(X \sim G)\) is the mixture of interventional distributions \(P(Y \mid do(X = x))\) weighted by \(G\).

Biological example: Heterogeneous treatment response, when the “same” treatment produces different effects across individuals, modelling \(do(X \sim G)\) captures the distribution of responses. For instance, a drug dose may vary around a target due to pharmacokinetic variability.

17.3.4 Dynamic Regimes

Dynamic regimes make the intervention depend on the current state (or history) of the system:

\[ do(X(t) = g(Y(t))) \]

where \(g\) is a function of the observed state \(Y(t)\) (or more generally, the history \(H_t\)). The intervention is not fixed in advance but adapts as the system evolves.

Feedback control as a dynamic regime: A controller that adjusts the input \(u(t)\) based on the current state \(x(t)\) implements a dynamic regime: \(u(t) = \pi(x(t))\).

Connection to optimal control theory: Optimal control seeks the policy \(\pi^*\) that maximises (or minimises) an objective. The control law \(\pi^*\) defines a dynamic regime.

Mathematical formulation:

\[ \dot{X} = f(X) + u(X, t) \]

where \(u\) is the control policy. The intervention \(u(X, t)\) depends on both state and time, making the system a closed-loop controlled system.

Biological example: Adaptive dosing protocols adjust drug dose based on measured biomarkers (e.g., therapeutic drug monitoring). The dose at time \(t\) depends on the patient’s current state rather than being fixed in advance.

17.3.5 Summary Table

Type Mathematical form When to use Biological example
Point \(X(t^*) = x^*\) or \(\dot{X} = f(X) + c \cdot \delta(t - t^*)\) Instantaneous, negligible duration Bolus drug injection
Sustained \(X(t) = x^*\) for \(t \in [t_1, t_2]\) Constant intervention over interval Continuous drug infusion
Stochastic \(do(X \sim G)\); \(P(Y \mid do(X \sim G)) = \int P(Y \mid do(X = x)) \, dG(x)\) Heterogeneous or uncertain intervention Heterogeneous treatment response
Dynamic regime \(do(X(t) = g(Y(t)))\); \(\dot{X} = f(X) + u(X, t)\) State-dependent, adaptive intervention Adaptive dosing protocols

17.4 Policy Evaluation in Deterministic Systems

17.4.1 Evaluating Treatment Policies

Question: What is the expected outcome under policy \(\pi\)?

\[ \mathbb{E}[Y^{do(\pi)}] = \mathbb{E}[Y \mid \text{policy } \pi] \]

Method: Simulate ODE under policy and compute expected outcomes.

17.4.2 Optimal Control

Question: What policy maximises expected outcome?

\[ \pi^* = \arg\max_\pi \mathbb{E}^{do(\pi)}[Y_T] \]

Methods: Dynamic programming, Pontryagin’s maximum principle.

17.4.3 Implementation: Policy Evaluation

We can evaluate policies by simulating the ODE under different treatment strategies:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))
@auto_using OrdinaryDiffEq CairoMakie

# Example: Treatment policy in a simple disease model
# State: X (disease severity), Action: A (treatment level 0 or 1)
# Policy π: Treat if severity > threshold

β = 0.2  # Disease progression rate
α = 0.3  # Treatment effectiveness

# Policy 1: No treatment (A = 0 always)
function disease_no_treatment!(du, u, p, t)
    X = u[1]
    du[1] = β * X  # Disease progresses
end

# Policy 2: Treat if X > 0.5 (threshold policy)
function disease_threshold!(du, u, p, t)
    X = u[1]
    A = X > 0.5 ? 1.0 : 0.0  # Policy: treat if severity > 0.5
    du[1] = β * X - α * A * X  # Treatment reduces severity
end

# Policy 3: Always treat (A = 1 always)
function disease_always_treat!(du, u, p, t)
    X = u[1]
    du[1] = β * X - α * X  # Always apply treatment
end

# Initial condition
u0 = [0.1]  # Initial severity = 0.1
tspan = (0.0, 20.0)

# Solve under each policy
prob1 = ODEProblem(disease_no_treatment!, u0, tspan)
sol1 = solve(prob1, Tsit5())

prob2 = ODEProblem(disease_threshold!, u0, tspan)
sol2 = solve(prob2, Tsit5())

prob3 = ODEProblem(disease_always_treat!, u0, tspan)
sol3 = solve(prob3, Tsit5())

# Evaluate expected outcomes (final severity)
outcome1 = sol1.u[end][1]
outcome2 = sol2.u[end][1]
outcome3 = sol3.u[end][1]

println("Expected final severity under each policy:")
println("  Policy 1 (no treatment): ", round(outcome1, digits=3))
println("  Policy 2 (threshold): ", round(outcome2, digits=3))
println("  Policy 3 (always treat): ", round(outcome3, digits=3))
println("\nBest policy: ", outcome3 < outcome2 < outcome1 ? "Policy 3 (always treat)" : 
        outcome2 < outcome1 ? "Policy 2 (threshold)" : "Policy 1 (no treatment)")
Expected final severity under each policy:
  Policy 1 (no treatment): 5.46
  Policy 2 (threshold): 0.502
  Policy 3 (always treat): 0.014

Best policy: Policy 3 (always treat)
Figure 17.2: Policy evaluation: comparing different treatment strategies

17.4.4 Optimal Control

For optimal control problems, we seek the policy that maximises (or minimises) an objective function. This typically requires specialised optimisation methods:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

# Note: Full optimal control implementation requires:
# - Optimization.jl for optimisation
# - Dynamic programming or Pontryagin's maximum principle
# - Cost function definition

# Example framework (conceptual):
# Objective: Minimise ∫[X(t)² + c*A(t)²] dt
#   - X(t)²: Penalty for high disease severity
#   - c*A(t)²: Cost of treatment (c is cost parameter)
# Constraints: 0 ≤ A(t) ≤ 1 (treatment bounded)

println("Optimal control framework:")
println("  Objective: Minimise ∫[X(t)² + c*A(t)²] dt")
println("  Subject to: dX/dt = β*X - α*A*X")
println("  Constraints: 0 ≤ A(t) ≤ 1")
println("\nMethods:")
println("  - Dynamic programming (discrete time)")
println("  - Pontryagin's maximum principle (continuous time)")
println("  - Numerical optimisation (Optimization.jl)")
println("\nFor full implementation, see Optimization.jl documentation")
Optimal control framework:
  Objective: Minimise ∫[X(t)² + c*A(t)²] dt
  Subject to: dX/dt = β*X - α*A*X
  Constraints: 0 ≤ A(t) ≤ 1

Methods:
  - Dynamic programming (discrete time)
  - Pontryagin's maximum principle (continuous time)
  - Numerical optimisation (Optimization.jl)

For full implementation, see Optimization.jl documentation

17.5 Stratum context

This chapter addresses Doing in the Dynamical stratum: how can we intervene/modify deterministic dynamic processes? Interventions in ODEs modify system dynamics, allowing us to evaluate policies and find optimal control strategies.

17.6 Key Takeaways

  1. Interventions modify ODEs: \(do(\cdot)\) operator changes system dynamics
  2. Policy evaluation: Simulate ODE under policy to evaluate outcomes
  3. Optimal control: Find policies that maximise expected outcomes
  4. Structural changes propagate: Interventions affect system evolution

17.7 Further Reading