diff --git a/CHANGELOG.md b/CHANGELOG.md
index 19373f8..b8dff90 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- **Time-varying load for `RegressionNode` (load-dependent aging, #37).**
+  `RegressionNode` now accepts a `schedule=` (a surpyval `StepSchedule`) as an
+  alternative to a fixed covariate vector: the load the component runs under
+  varies over its life, and reliability is the exact survival along that
+  covariate path (`model.sf_tvc`). Conditioning on `age` gives the go-forward
+  reliability from the component's current life under the schedule — the
+  digital-twin / load-dependent-aging node — with no change to the
+  condition-based layer, since `sf_tvc(age+x)/sf_tvc(age)` is exactly surpyval's
+  `sf_tvc(..., given=age)`. Works for accelerated-failure-time and
+  proportional-/additive-hazards families (not proportional-odds). The schedule
+  persists through serialisation. (Uses surpyval's `sf_tvc`, which is why the
+  minimum surpyval is now 0.16 — see below.)
 - **`LoadSharingModel`: load-sharing dynamic node (dependent failure).** A
   sibling to `StandbyModel` where *n* coupled units share a total load and the
   survivors carry more (and so age faster) as siblings fail — the group works
@@ -35,6 +47,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   (`structural_importance`, being probability-free, is unaffected).
   Alpha-factor is a planned extension. (#44)
 
+### Changed
+- Require **surpyval >= 0.16** (was >= 0.15): the time-varying-load
+  `RegressionNode` schedule mode is built on surpyval's `sf_tvc`, added in 0.16.
+
+### Documentation
+- Concepts gains theory sections for covariate/time-varying reliability, load
+  sharing, and common-cause failures; the tutorial gains a worked
+  dependent-failures step exercising all three. (#37, #38, #44)
+
 ## [0.7.0] - 2026-07-20
 
 The **Maintenance & Covariates** milestone: price imperfect-repair (generalized
diff --git a/docs/concepts.md b/docs/concepts.md
index 15db618..d7b188f 100644
--- a/docs/concepts.md
+++ b/docs/concepts.md
@@ -116,6 +116,131 @@ is a cheap, exact re-evaluation — the heterogeneous generalisation of the
 age, to *per-component* ages. Only lifetime (time-varying) distributions age; a
 fixed-probability component's reliability does not depend on `Xᵢ`.
 
+## Covariate-dependent and time-varying reliability
+
+A component's reliability often depends on the *conditions it runs under*, not
+only on elapsed time. If you have fitted a **regression** survival model in
+surpyval — accelerated-failure-time (AFT), proportional-hazards (Cox, PH),
+proportional-odds (PO), … — a [`RegressionNode`][repyability.RegressionNode]
+uses it as an ordinary RBD node.
+
+**Fixed covariates.** Pin the component's operating point `Z` (temperature,
+load, duty cycle) and its reliability is the model's survival there:
+
+```
+R(x) = model.sf(x, Z)
+```
+
+That is a single univariate curve, so the node takes part in system reliability,
+importance, MTTF and the condition-based (`age`) layer with no special handling
+— a hotter-running unit is just a node with different covariates. This is
+family-agnostic: AFT, PH and PO all expose `sf(x, Z)`.
+
+**A time-varying covariate path.** When the load itself changes over the
+component's life, the reliability is no longer the survival at one covariate
+value but the survival *along the whole path* `Z(t)`. For a piecewise-constant
+path (a surpyval `StepSchedule`) this is
+
+```
+R(x) = model.sf_tvc(x, schedule)
+```
+
+the probability of surviving each segment in turn under its own covariate.
+Conditioning on an `age` needs no special case, because
+
+```
+R(x | age) = sf_tvc(age + x) / sf_tvc(age) = sf_tvc(x, schedule, given=age)
+```
+
+is exactly the go-forward survival from the component's current life under the
+schedule — the load-dependent-aging / "digital twin" node. Whether a family
+composes along a path is a property of the model: **AFT** (the path rescales the
+clock) and **proportional-/additive-hazards** (the path accumulates hazard) do;
+**proportional-odds** does not, and is refused in schedule mode. The
+fixed-covariate node is the special case of a constant path.
+
+## Dependent failures: load sharing
+
+Redundant units that *share a load* do not fail independently. While all are up
+each carries its share; when one fails the survivors pick up the slack, run
+harder, and age faster — so the failures are positively correlated, and treating
+them as `n` independent parallel nodes over-counts the redundancy.
+
+A [`LoadSharingModel`][repyability.LoadSharingModel] captures the coupling as a
+single node. Each unit is a fitted AFT model with **load as its covariate**, so a
+unit running under load `ℓ` ages on its baseline clock at an acceleration factor
+`φ(ℓ)` (the AFT time-scaling). With `s` survivors sharing a total load `L`, each
+carries `L / s` and ages at `φ(L / s)`; as siblings fail, `s` falls, `L / s`
+rises, and the survivors' clocks speed up. The group works while at least `k` of
+the `n` units survive. This is the **cumulative-exposure** model: a unit's
+*virtual age* is the integral of `φ(load(t))` over real time, and it fails when
+that virtual age reaches its baseline failure age.
+
+Two regimes:
+
+- **Closed form.** Identical units with an **Exponential** baseline give a group
+  lifetime that is a sum of exponential stages — each stage the time for the
+  next unit to fail at the current shared load — i.e. a **hypoexponential**
+  distribution, evaluated exactly with no simulation (`is_simulated == False`).
+- **Simulation.** Otherwise the survival curve is a Kaplan–Meier fit to
+  lifetimes drawn from the cumulative-exposure event loop (seeded;
+  `is_simulated == True`).
+
+As a check, with no load effect (`φ ≡ 1`) the survivors do not accelerate and the
+group reduces *exactly* to the ordinary k-out-of-n parallel result. Load sharing
+is the "self-loading" sibling of the condition-based layer: there the load is
+streamed in from sensors, here it is computed from the group's own survivors.
+
+## Common-cause failures
+
+Redundancy only buys reliability if the redundant units fail for *independent*
+reasons. In practice they often share a root cause — a common manufacturing
+batch, a shared power supply, one miscalibration applied to every unit — and a
+single event takes them all down together. Because the exact engine assumes
+independence, it **over-estimates** a redundant group; a common-cause model
+injects the shared coupling. (This is the mirror image of load sharing: there the
+coupling is mechanical load transfer, here it is a shared shock.)
+
+A [`CCFGroup`][repyability.CCFGroup] declares the coupled (symmetric) members and
+the model, passed via `ccf_groups`. Two models:
+
+- [`BetaFactor(beta)`][repyability.BetaFactor] — a fraction `β` of each unit's
+  failure probability `Q` comes from a cause shared across the **whole** group
+  (which fails every member at once); the remaining `(1 − β) Q` is independent.
+  The workhorse of probabilistic-risk assessment.
+- [`MGL(beta, gamma, ...)`][repyability.MGL] — the **Multiple Greek Letter**
+  model, which also resolves *partial* common causes (a cause failing some but
+  not all of the group) through a cascade of conditional probabilities:
+  `β = P(shared by ≥ 2 | failed)`, `γ = P(≥ 3 | ≥ 2)`, and so on. The probability
+  that a cause fails a *specific* set of `k` of the `m` members is the standard
+  basic-event probability
+
+```
+Q_k = [ 1 / C(m−1, k−1) ] · (ρ₁ ρ₂ ⋯ ρ_k) · (1 − ρ_{k+1}) · Q
+```
+
+  with `ρ₁ = 1, ρ₂ = β, ρ₃ = γ, …, ρ_{m+1} = 0`; these partition each unit's `Q`
+  exactly. A group of `m` members takes `m − 1` letters, and `MGL(β)` on two
+  members is exactly `BetaFactor(β)`.
+
+**The evaluation is exact, not a correction factor.** Each model's `decompose`
+splits the group's failure into an independent part plus a set of
+**mutually-exclusive shocks** (each a subset of members failing together). The
+system reliability is then computed by conditioning on the shock outcome of every
+group — in each branch the shocked members are down and the rest fail only
+independently, so it is an ordinary independent system-reliability evaluation —
+and blending the branches by their probabilities. Hence `β = 0` reproduces the
+independent result exactly, and `β = 1` makes a redundant group no better than a
+single unit.
+
+Common cause is currently reflected in `sf()` / `ff()` (and quantities derived
+from them) and persists through serialisation. Groups must be symmetric
+(identical member models) and disjoint. The probability-dependent
+importance/sensitivity and the condition-based methods do not yet account for it
+and raise a clear error on a CCF RBD; `structural_importance`, being
+probability-free, is unaffected. **Alpha-factor**, a data-estimable
+reparameterisation of the same multiplicities, is a planned extension.
+
 ## Scope
 
 - **Fitting failure data to distributions lives in
diff --git a/docs/guide.md b/docs/guide.md
index 183c5cc..7850f7d 100644
--- a/docs/guide.md
+++ b/docs/guide.md
@@ -227,6 +227,22 @@ rbd.remaining_life(0.9, {"m": NodeState(age=40)})
   meaning (operating time survived): `R(x | age) = sf(age + x, Z) / sf(age, Z)`,
   which holds for every regression family. A hotter-running motor is simply a
   node with different covariates.
+- **Time-varying load (a `schedule`).** If the load *changes over the
+  component's life*, pass a surpyval `StepSchedule` instead of a fixed vector —
+  `RegressionNode(model, schedule=...)`. Reliability is then the exact survival
+  along that covariate path (`model.sf_tvc`), and conditioning on `age` gives
+  the go-forward reliability from the component's current life under the
+  schedule — the load-dependent-aging / digital-twin model of issue #37. This
+  works for accelerated-failure-time and proportional-/additive-hazards models
+  (not proportional-odds), and needs a surpyval build that provides `sf_tvc`.
+
+```python
+from surpyval import StepSchedule
+# benign for 500 h, then a harsher load
+sched = StepSchedule.from_changepoints([0, 500], [[0.0], [1.0]])
+motor = RegressionNode(aft, schedule=sched)
+rbd.sf_given_state(x, {"m": NodeState(age=520)})  # forward reliability, now in the harsh phase
+```
 - The node serialises with the RBD — the fitted regression model round-trips
   through `surpyval.from_dict`.
 - `mean` and `random` (simulation-based MTTF) need a proper parametric lifetime,
diff --git a/docs/index.md b/docs/index.md
index b30a792..2b6aaf5 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -109,5 +109,6 @@ rbd.mean_availability()
 - **[User guide](guide.md)** — building RBDs, reliability, importance measures,
   forcing nodes, condition-based evaluation, and availability results.
 - **[Concepts](concepts.md)** — the theory: path/cut sets, choosing an
-  importance measure, and how conditioning works.
+  importance measure, how conditioning works, and the dependent-failure models
+  (load sharing and common cause) and covariate/time-varying reliability.
 - **[API reference](api.md)** — every public class and method.
diff --git a/docs/tutorial.md b/docs/tutorial.md
index 088179e..3b53737 100644
--- a/docs/tutorial.md
+++ b/docs/tutorial.md
@@ -2,9 +2,9 @@
 
 This walkthrough takes a small but realistic system and answers the questions a
 reliability engineer actually asks of it — *how reliable is it, when should we
-service it, what is the weak link, and what does its remaining life look like
-once it is in the field.* Every capability in RePyability shows up here in the
-order you would reach for it.
+service it, what is the weak link, what does its remaining life look like once it
+is in the field, and where is its redundancy weaker than it looks.* Every
+capability in RePyability shows up here in the order you would reach for it.
 
 We will model a **pumping skid**:
 
@@ -184,7 +184,92 @@ Stream a new `state` each time fresh readings arrive; every update is a cheap,
 exact re-evaluation. (Only lifetime distributions age; standby/composite nodes
 are out of scope for the condition-based methods in this release.)
 
-## 8. If the skid is repairable: availability
+## 8. When redundancy is optimistic: dependent failures
+
+Steps 1–7 assumed the five components fail **independently** and each runs at a
+**fixed** load. Both assumptions flatter a redundant system. RePyability lets you
+relax each one where it matters — and because the system quantity stays exact,
+you see precisely what the optimism was worth. The theory behind all three is in
+[Concepts](concepts.md).
+
+### The pumps share the load
+
+The parallel-pump model assumed the surviving pump is unaffected when its sibling
+drops. In reality it then carries the *whole* duty, runs harder, and ages faster.
+Refit the pump with **load as a covariate** and model the pair as a load-sharing
+group — each pump carries half of a total load of `2.0` while both run, and the
+group needs at least one:
+
+```python
+import surpyval as surv
+from repyability import LoadSharingModel
+
+pump = surv.ExponentialAFT.fit(pump_hours, Z=pump_loads)   # fitted in surpyval
+pumps = LoadSharingModel([pump, pump], load=2.0, k=1)       # 2 units, need >= 1
+
+pumps.sf(4000)     # -> 0.868   the pair, sharing the load
+```
+
+Modelled *independently* — the same two pumps each pinned at the half-load — the
+redundant pair would read `0.920`. The load transfer has quietly eaten a third of
+the pair's redundancy margin, and more as it wears: `0.641` vs `0.764` at 8000 h.
+Identical Exponential-baseline units like these get the **exact** hypoexponential
+group lifetime, so `pumps.is_simulated` is `False` — no simulation involved.
+
+### A common cause on the filters
+
+The two filters are the same part from the same shelf, so one bad batch — or one
+contamination event upstream — can blind both at once. That is a **common-cause
+failure**, and no amount of *structural* redundancy defends against it. Attach a
+beta-factor group (here 8% of a filter's failures are shared) to the otherwise
+unchanged skid:
+
+```python
+from repyability import CCFGroup, BetaFactor
+
+rbd_ccf = NonRepairableRBD(
+    edges, reliabilities,
+    ccf_groups=[CCFGroup(["filterA", "filterB"], BetaFactor(0.08))],
+)
+
+rbd_ccf.sf(4000)   # -> 0.9023   vs 0.9091 with independent filters
+```
+
+Modest at the skid level *here*, because the filters are not the weak link — but
+the coupling is exact, and on a system that leans on its redundant pair it is the
+difference between a design that meets its target and one that only appears to.
+`BetaFactor(0)` recovers the independent number; `MGL(β, γ, …)` handles a cause
+that fails *some* but not all of a larger group.
+
+### A duty that ramps up
+
+Finally, drop the fixed-load assumption. Say the skid is commissioned gently and
+then, at 3000 h, pushed to a harsher continuous duty. Fit a regression model with
+the duty as a covariate and give the node a **schedule** instead of a single
+covariate vector:
+
+```python
+from surpyval import StepSchedule
+from repyability import RegressionNode
+
+duty_unit = surv.WeibullAFT.fit(run_hours, Z=duty_history)   # fitted in surpyval
+# benign until 3000 h, then a harsher duty for the rest of life
+duty = StepSchedule.from_changepoints([0, 3000], [[0.0], [1.0]])
+node = RegressionNode(duty_unit, schedule=duty)
+
+node.sf(4000)      # -> 0.822   just after the step up
+node.sf(6000)      # -> 0.603   the harsher duty has now done real damage
+```
+
+Held at the benign duty the same unit would read `0.861` and `0.736`. Reliability
+is the exact survival *along* the duty path, and conditioning on the component's
+age (`NodeState(age=...)`) gives its go-forward reliability from wherever it sits
+on that path — the digital twin of step 7, now with a load history.
+
+Each of these is an ordinary node: it serialises with the RBD and takes part in
+the same exact system computation as everything above.
+
+## 9. If the skid is repairable: availability
 
 If a failed component is repaired rather than replaced, give each one a
 *repairability* (time-to-repair) distribution and ask about availability instead
@@ -212,7 +297,7 @@ Long-run availability is exact and needs no simulation; the time-resolved
 availability curve and its criticality measures come from a seeded
 discrete-event simulation.
 
-## 9. Persist the model
+## 10. Persist the model
 
 The structure and its fitted models round-trip to plain JSON, so you can build
 the RBD once and reload it wherever the analysis runs (a dashboard, a scheduled
@@ -229,8 +314,9 @@ transient input you supply at evaluation time.
 ## Where to next
 
 - **[Concepts](concepts.md)** — the theory behind these numbers: path/cut sets,
-  the full importance-measure family and when to use each, and how conditioning
-  works.
+  the full importance-measure family and when to use each, how conditioning
+  works, and the dependent-failure models from step 8 (load sharing, common
+  cause, and covariate/time-varying reliability).
 - **[User guide](guide.md)** — the reference for every method, argument, and
   return contract.
 - **[API reference](api.md)** — the generated signatures and docstrings.
diff --git a/pyproject.toml b/pyproject.toml
index 8fa2f04..f21db76 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -35,7 +35,7 @@ dependencies = [
     "numpy>=2.0,<3",
     "scipy>=1.7",
     "networkx>=3.0",
-    "surpyval>=0.15,<0.16",
+    "surpyval>=0.16,<0.17",
     "tqdm>=4.64",
 ]
 
diff --git a/repyability/rbd/regression_node.py b/repyability/rbd/regression_node.py
index b646894..f5b9ac9 100644
--- a/repyability/rbd/regression_node.py
+++ b/repyability/rbd/regression_node.py
@@ -2,42 +2,53 @@
 
 A :class:`RegressionNode` wraps a fitted surpyval **regression** model (an
 accelerated-failure-time, proportional-hazards, proportional-odds, ... model)
-together with a fixed **covariate vector** ``Z`` -- the operating conditions of
-this component in the system. Its reliability is simply the model's survival at
-those covariates::
-
-    R(x) = model.sf(x, Z)
-
-so a regression node is an ordinary univariate node: it takes part in system
-reliability, importance, MTTF and the condition-based (``age``) layer exactly
-like any other node, with no special handling. Conditioning on age keeps its
-usual meaning -- operating time survived -- because at a fixed covariate the
-component has a plain univariate lifetime and
-``R(x | age) = sf(age + x, Z) / sf(age, Z)`` holds for every regression family.
-
-RePyability *consumes* the fitted model; do the regression fit in surpyval.
-See issue #37.
+together with the covariate history of *this* component in the system, and
+exposes its survival as an ordinary univariate node. The covariates can be:
+
+* a **fixed vector** ``Z`` (constant operating conditions) -- reliability is
+  ``R(x) = model.sf(x, Z)``; or
+* a **time-varying schedule** ``Z(t)`` (a surpyval ``StepSchedule``: the load
+  the component runs under changes over its life) -- reliability is
+  ``R(x) = model.sf_tvc(x, schedule)``, the exact survival along that
+  covariate path (accelerated-failure-time / proportional- / additive-hazards;
+  not proportional-odds). This is the load-dependent-aging / digital-twin node
+  of issue #37: as-new survival integrates the whole load path, and
+  conditioning on ``age`` gives the go-forward reliability from the component's
+  current life, since ``R(x | age) = sf_tvc(age + x) / sf_tvc(age)`` is exactly
+  surpyval's ``sf_tvc(..., given=age)``.
+
+Either way it is an ordinary univariate node -- it takes part in system
+reliability, importance, MTTF and the condition-based (``age``) layer with no
+special handling. RePyability *consumes* the fitted model; do the regression
+fit in surpyval.
 """
 
-from typing import Any
+from typing import Any, Optional
 
 import numpy as np
 from numpy.typing import ArrayLike
 
 
 class RegressionNode:
-    """An RBD node backed by a fitted surpyval regression model at fixed
-    covariates.
+    """An RBD node backed by a fitted surpyval regression model.
+
+    Provide exactly one of ``covariates`` (a fixed operating point) or
+    ``schedule`` (a time-varying covariate path).
 
     Parameters
     ----------
     model : surpyval regression model
-        A fitted regression model whose ``sf(x, Z)`` gives survival at a
-        covariate matrix ``Z`` (e.g. ``surpyval.WeibullAFT.fit(...)``,
-        ``surpyval.CoxPH.fit(...)``).
-    covariates : array_like
-        The component's covariate vector ``Z`` (the operating conditions),
-        matching the covariates the model was fitted with.
+        A fitted regression model, e.g. ``surpyval.WeibullAFT.fit(...)`` or
+        ``surpyval.CoxPH.fit(...)``. Fixed covariates use its ``sf(x, Z)``; a
+        schedule uses its ``sf_tvc(x, schedule)`` (needs a surpyval that
+        provides it, and a family other than proportional-odds).
+    covariates : array_like, optional
+        The component's fixed covariate vector ``Z`` (its operating
+        conditions), matching the covariates the model was fitted with.
+    schedule : surpyval StepSchedule, optional
+        A piecewise-constant covariate path ``Z(t)`` -- the load the component
+        runs under over its life (build with
+        ``surpyval.StepSchedule.from_changepoints`` / ``from_intervals``).
 
     Examples
     --------
@@ -53,41 +64,65 @@ class RegressionNode:
     True
     """
 
-    def __init__(self, model: Any, covariates: ArrayLike):
+    def __init__(
+        self,
+        model: Any,
+        covariates: Optional[ArrayLike] = None,
+        schedule: Any = None,
+    ):
+        if (covariates is None) == (schedule is None):
+            raise ValueError(
+                "RegressionNode requires exactly one of `covariates` (a fixed "
+                "covariate vector) or `schedule` (a time-varying "
+                "StepSchedule)."
+            )
         self.model = model
-        self.covariates = np.atleast_1d(np.asarray(covariates, dtype=float))
-        # Probe the regression interface so a misuse (a non-regression model,
-        # or covariates of the wrong width) fails clearly at construction.
+        self.covariates = (
+            None
+            if covariates is None
+            else np.atleast_1d(np.asarray(covariates, dtype=float))
+        )
+        self.schedule = schedule
+        # Probe the survival interface so a misuse fails clearly at
+        # construction (wrong covariate width, an sf_tvc-less surpyval, or a
+        # proportional-odds model in schedule mode).
         try:
             probe = self._sf_at(np.array([1.0]))
             if not np.all(np.isfinite(probe)):
-                raise ValueError("sf(x, Z) returned non-finite values")
+                raise ValueError("survival returned non-finite values")
         except Exception as e:
             raise ValueError(
-                "RegressionNode requires a fitted surpyval regression model "
-                "whose sf(x, Z) accepts a covariate matrix, and covariates "
-                f"matching the model's fitted width. Probing sf failed: "
+                "RegressionNode requires a fitted surpyval regression model. "
+                "In fixed-covariate mode its sf(x, Z) must accept a covariate "
+                "matrix of the fitted width; in schedule mode the model must "
+                "support sf_tvc(x, schedule) (accelerated-failure-time / "
+                "proportional- or additive-hazards on a recent surpyval, not "
+                f"proportional-odds). Probing survival failed: "
                 f"{type(e).__name__}: {e}."
             ) from e
         # Cached (t, sf(t)) grid for mean()/random() (built lazily).
         self._grid: Any = None
 
     def _Z(self, n: int) -> np.ndarray:
-        """The covariate vector broadcast to ``n`` rows for ``sf(x, Z)``."""
+        """Fixed covariate vector broadcast to ``n`` rows for ``sf(x, Z)``."""
+        assert self.covariates is not None  # fixed-covariate mode only
         return np.repeat(self.covariates[np.newaxis, :], n, axis=0)
 
     def _sf_at(self, x: np.ndarray) -> np.ndarray:
         x = np.atleast_1d(np.asarray(x, dtype=float))
+        if self.schedule is not None:
+            return np.asarray(self.model.sf_tvc(x, self.schedule), dtype=float)
         return np.asarray(self.model.sf(x, self._Z(len(x))), dtype=float)
 
     # -- Node reliability interface ---------------------------------------
 
     def sf(self, x: ArrayLike) -> np.ndarray:
-        """Reliability at the stored covariates: ``model.sf(x, Z)``."""
+        """Reliability at the stored covariates (``model.sf(x, Z)``) or along
+        the covariate schedule (``model.sf_tvc(x, schedule)``)."""
         return self._sf_at(np.atleast_1d(np.asarray(x, dtype=float)))
 
     def ff(self, x: ArrayLike) -> np.ndarray:
-        """Unreliability at the stored covariates: ``1 - sf(x)``."""
+        """Unreliability: ``1 - sf(x)``."""
         return 1.0 - self.sf(x)
 
     def _survival_grid(self):
@@ -124,7 +159,7 @@ def _survival_grid(self):
         return self._grid
 
     def mean(self) -> float:
-        """Mean time to failure at the stored covariates.
+        """Mean time to failure at the stored covariates / along the schedule.
 
         For a non-negative lifetime ``E[T] = integral of R(t)``, integrated
         numerically over the survival curve.
@@ -133,7 +168,7 @@ def mean(self) -> float:
         return float(np.trapezoid(s, t))
 
     def random(self, size: int) -> np.ndarray:
-        """Draw ``size`` failure times at the stored covariates.
+        """Draw ``size`` failure times at the stored covariates / schedule.
 
         Inverse-transform sampling on the survival curve. Uses numpy's global
         RNG, so wrap the call in
@@ -146,13 +181,26 @@ def random(self, size: int) -> np.ndarray:
 
     # -- Serialisation ----------------------------------------------------
 
+    @staticmethod
+    def _schedule_to_dict(schedule: Any) -> dict:
+        if getattr(schedule, "period", None) is not None:
+            raise NotImplementedError(
+                "Serialising a cyclic StepSchedule is not supported yet; use "
+                "a change-point / interval schedule."
+            )
+        times = [float(e) for e in schedule.edges if np.isfinite(e)]
+        return {"times": times, "values": np.asarray(schedule.Z).tolist()}
+
     def to_dict(self) -> dict:
-        """Serialise to a JSON-friendly dict (the fitted model + covariates).
-        See :meth:`from_dict`."""
-        return {
-            "model": self.model.to_dict(),
-            "covariates": [float(v) for v in self.covariates],
-        }
+        """Serialise to a JSON-friendly dict (the fitted model + covariates or
+        schedule). See :meth:`from_dict`."""
+        out: dict = {"model": self.model.to_dict()}
+        if self.schedule is not None:
+            out["schedule"] = self._schedule_to_dict(self.schedule)
+        else:
+            assert self.covariates is not None
+            out["covariates"] = [float(v) for v in self.covariates]
+        return out
 
     @classmethod
     def from_dict(cls, d: dict) -> "RegressionNode":
@@ -160,10 +208,21 @@ def from_dict(cls, d: dict) -> "RegressionNode":
         through ``surpyval.from_dict``)."""
         import surpyval
 
-        return cls(surpyval.from_dict(d["model"]), d["covariates"])
+        model = surpyval.from_dict(d["model"])
+        if "schedule" in d:
+            from surpyval.univariate.regression import StepSchedule
+
+            sd = d["schedule"]
+            schedule = StepSchedule.from_changepoints(
+                sd["times"], sd["values"]
+            )
+            return cls(model, schedule=schedule)
+        return cls(model, covariates=d["covariates"])
 
     def __repr__(self) -> str:
-        return (
-            f"RegressionNode({type(self.model).__name__}, "
-            f"covariates={list(self.covariates)})"
-        )
+        if self.schedule is not None:
+            which = f"schedule={self.schedule!r}"
+        else:
+            cov = list(self.covariates)  # type: ignore[arg-type]
+            which = f"covariates={cov}"
+        return f"RegressionNode({type(self.model).__name__}, {which})"
diff --git a/repyability/tests/test_regression_node.py b/repyability/tests/test_regression_node.py
index d02e209..77dbc7d 100644
--- a/repyability/tests/test_regression_node.py
+++ b/repyability/tests/test_regression_node.py
@@ -16,6 +16,20 @@
 
 from repyability import NodeState, NonRepairableRBD, RegressionNode
 
+# The time-varying-covariate (schedule) mode needs a surpyval build exposing
+# sf_tvc / StepSchedule; skip those tests on older builds so CI stays green.
+try:
+    from surpyval.univariate.regression import StepSchedule
+
+    _HAS_TVC = True
+except ImportError:  # pragma: no cover
+    StepSchedule = None
+    _HAS_TVC = False
+
+needs_tvc = pytest.mark.skipif(
+    not _HAS_TVC, reason="surpyval build lacks sf_tvc / StepSchedule"
+)
+
 
 @pytest.fixture(scope="module")
 def data():
@@ -182,3 +196,97 @@ def test_non_regression_model_rejected():
         RegressionNode(
             surv.Weibull.from_params([100.0, 2.0]), covariates=[0.1]
         )
+
+
+def test_requires_exactly_one_of_covariates_or_schedule(models):
+    with pytest.raises(ValueError, match="exactly one"):
+        RegressionNode(models["aft"])  # neither
+
+
+# -- time-varying covariate (schedule) mode, issue #37 ---------------------
+
+
+@needs_tvc
+@pytest.mark.parametrize("family", ["aft", "ph"])
+def test_schedule_sf_matches_sf_tvc(models, family):
+    sched = StepSchedule.from_changepoints([0, 50], [[0.0], [0.8]])
+    node = RegressionNode(models[family], schedule=sched)
+    xt = np.array([30.0, 80.0, 150.0])
+    assert np.allclose(node.sf(xt), np.ravel(models[family].sf_tvc(xt, sched)))
+
+
+@needs_tvc
+def test_schedule_condition_based_is_given(models):
+    # The condition-based `age` path must equal surpyval's sf_tvc(given=age):
+    # forward reliability from the component's current life under the schedule.
+    sched = StepSchedule.from_changepoints([0, 50], [[0.0], [0.8]])
+    node = RegressionNode(models["aft"], schedule=sched)
+    rbd = NonRepairableRBD([("s", "c"), ("c", "t")], {"c": node})
+    age = 60.0
+    for x in (20.0, 40.0):
+        got = float(rbd.sf_given_state(x, {"c": NodeState(age=age)}))
+        want = float(
+            np.ravel(
+                models["aft"].sf_tvc(np.array([age + x]), sched, given=age)
+            )[0]
+        )
+        assert got == pytest.approx(want)
+
+
+@needs_tvc
+def test_schedule_more_load_lowers_reliability(models):
+    # A schedule ramping to a higher load is less reliable than a benign one.
+    low = RegressionNode(
+        models["aft"],
+        schedule=StepSchedule.from_changepoints([0, 50], [[0.0], [0.0]]),
+    )
+    high = RegressionNode(
+        models["aft"],
+        schedule=StepSchedule.from_changepoints([0, 50], [[0.0], [1.0]]),
+    )
+    xt = np.array([120.0])
+    assert high.sf(xt)[0] < low.sf(xt)[0]
+
+
+@needs_tvc
+def test_po_schedule_rejected(models):
+    # Proportional odds does not compose along a time-varying covariate.
+    with pytest.raises(ValueError, match="proportional-odds|sf_tvc"):
+        RegressionNode(
+            models["po"],
+            schedule=StepSchedule.from_changepoints([0, 50], [[0.0], [0.8]]),
+        )
+
+
+@needs_tvc
+def test_schedule_serialisation_roundtrip(models):
+    sched = StepSchedule.from_changepoints([0, 50, 120], [[0.0], [0.8], [0.3]])
+    node = RegressionNode(models["aft"], schedule=sched)
+    node2 = RegressionNode.from_dict(json.loads(json.dumps(node.to_dict())))
+    xt = np.array([30.0, 90.0, 160.0])
+    assert np.allclose(node.sf(xt), node2.sf(xt))
+
+
+@needs_tvc
+def test_rbd_with_schedule_node_json_roundtrip(models):
+    sched = StepSchedule.from_changepoints([0, 50], [[0.0], [0.8]])
+    rbd = NonRepairableRBD(
+        [("s", "c"), ("c", "t")],
+        {"c": RegressionNode(models["aft"], schedule=sched)},
+    )
+    restored = NonRepairableRBD.from_json(rbd.to_json())
+    assert np.isclose(float(rbd.sf(80.0)), float(restored.sf(80.0)))
+
+
+@needs_tvc
+def test_schedule_mean_and_random(models):
+    node = RegressionNode(
+        models["aft"],
+        schedule=StepSchedule.from_changepoints([0, 50], [[0.0], [0.5]]),
+    )
+    assert node.mean() > 0.0
+    np.random.seed(0)
+    a = node.random(2000)
+    np.random.seed(0)
+    b = node.random(2000)
+    assert np.allclose(a, b) and (a > 0).all()