Reaction networks and differential/algebraic equations

These blocks bind reusable definitions, not trajectories. Assignment is the canonical spelling; REACTIONS name … END and EQUATIONS name … END are equivalent alternatives. Reading a definition returns its schema and model-local name without running integration. Solve calls are ordinary reactive formulas and never write back to the declared species or states. Begin either block form with a species, state, or parameter declaration.

input transfer_rate = 0.1
 
conversion = REACTIONS
  SPECIES A INITIAL 10
  SPECIES B INITIAL 0
  PARAMETER k = transfer_rate
  REACTION transfer: A -> B AT k*A
END
 
trajectory = REACTION_SOLVE(conversion, 0, 20, {samples: 101})
columns = REACTION_SPECIES(conversion)

The trajectory has time in column 1, followed by species in declaration order. columns is [["A", "B"]]. A changed transfer_rate invalidates the trajectory; parameters and initial values are captured once for each solve. Parameters can use earlier parameters and previously declared numeric model bindings. Initial values can use earlier parameters or model bindings. Rates can use species, parameters, TIME, arithmetic and the supported pure numeric builtin subset. Capture outer values with PARAMETER; do not refer to them directly in rates.

A rate is the complete authored flux. 2 A -> B AT k*A means A' = -2*k*A, B' = k*A; it does not silently square A. Repeated species terms are combined with exact integer stoichiometry. Use 0 -> A for a source, A -> 0 for a sink. Each reaction rate is evaluated once per right-hand-side evaluation, in source order, and reused across species. Signed rates are allowed. There is no concentration clamp or automatic conservation correction.

For a reversible reaction, give two complete fluxes, forward then reverse:

REACTION binding: enzyme + substrate <-> complex AT k_on*enzyme*substrate, k_off*complex

AT MASS_ACTION(k) explicitly multiplies k by each reactant raised to its stoichiometric coefficient. 2 A -> B AT MASS_ACTION(k) means k*A^2; A + A -> B has the same exponent. With <->, each direction uses its own reactant side: A <-> B AT MASS_ACTION(k_on), MASS_ACTION(k_off) means k_on*A forward and k_off*B reverse. A source with no reactants uses k alone. This is deterministic concentration/amount kinetics, with no stochastic combinatorial factor. MASS_ACTION is contextual to the complete AT rate; it is not an ordinary callable or an implicit default for an authored flux.

Use block-local LET for a shared rate, force, or other intermediate:

LET saturation = maximum_rate * substrate / (affinity + substrate)
REACTION conversion: substrate -> product AT saturation

Declare LET values before all reactions/equations. They can read states, parameters, TIME, and earlier LET values. Forward references between LET values, cycles, and shadowing another local are rejected. Each is evaluated and checked for a finite numeric result once per solver callback, even when unused. PARAMETER values instead remain fixed throughout one solve. These are private model intermediates; output columns still contain declared species/states only.

charging = EQUATIONS
  PARAMETER capacitance = 0.1
  STATE voltage INITIAL 0
  capacitance * DERIVATIVE(voltage) = 5 - voltage
END
 
trajectory = EQUATION_SOLVE(charging, 0, 1, {samples: 101})
columns = EQUATION_STATES(charging)

Grid isolates diagonal derivative coefficients that are constant during a solve and calls ODE_SOLVE. For a constant mass matrix and algebraic constraints it calls DAE_SOLVE:

circuit = EQUATIONS
  STATE voltage INITIAL 0
  ALGEBRAIC current INITIAL 0
  0.1 * DERIVATIVE(voltage) = current
  0 = current - (5 - voltage)
END
 
trajectory = EQUATION_SOLVE(circuit, 0, 1, {samples: 101})

Write exactly one equation per variable, in declaration order. A STATE row has derivatives; an ALGEBRAIC row has none. Algebraic variables cannot appear inside DERIVATIVE. Their INITIAL values are guesses: the native DAE solver projects them to a consistent initial state. The differential mass block must be nonsingular. The supported DAE profile is semi-explicit, constant-mass, index 1; state-dependent mass, nonlinear derivatives, higher-index constraints, index reduction and purely algebraic systems are unsupported. Coefficients of derivatives may use finite constants and captured parameters, such as mass, 2*mass, or 1/capacitance. Their arithmetic and finite checks run once before integration. They cannot depend on states, algebraic variables, TIME, or per-callback LET values. A parameter update that makes a differential row zero or its mass block singular fails; it never reclassifies a STATE as ALGEBRAIC.

The solve, final, and report functions accept (definition, start, end[, options]). Use a direct same-module definition name or an immutable ordinary alias. Arbitrary records, strings, computed selections and imported definition references are unsupported in this initial profile. Definitions cannot be reassigned or written by rules. Assignment decorators and cell/qualified targets are unsupported for blocks. Solve calls may appear in same-module callables and local expressions. Local bindings must not shadow the selected definition or its captured dependencies; ambiguous captures are rejected at compilation.

The literal options record can also contain named scenario overrides:

baseline = REACTION_FINAL(conversion, 0, 20)
faster = REACTION_FINAL(conversion, 0, 20, {parameters: {k: 0.2}, initial: {A: 5}})

parameters overrides declared PARAMETER values; initial overrides declared species/state/algebraic initial values. Keys match local names without case sensitivity. Unknown or duplicate names and nonliteral override records are rejected. Values use pure numeric expressions in the caller's scope, so a local callable can pass a scenario value. Overrides do not bind the model's local names in caller scope. Default parameters/initials that reference earlier parameters see their overridden values. Each solve is independent and does not modify the definition, another solve, or the reactive inputs. The same finite and reaction-initial nonnegativity checks apply to overrides.

REACTION_FINAL and EQUATION_FINAL return the terminal scalar or a single state row. REACTION_REPORT and EQUATION_REPORT return the selected native solver's diagnostics. ODE reports support {format: "record"} with named work counters and terminal state; DAE reports retain their numeric row with initial projection and algebraic-residual diagnostics. Native report layouts differ. Each explicit solve/final/report call requests its own result; requesting three outputs does not promise one shared integration. Ordinary formula caching still applies.

Remaining options use the existing native solver's literal record, defaults, validation, budgets and cancellation. ODE methods are auto, rk45, or stiff; DAEs use the native implicit method. Unknown keys, invalid samples/tolerances and exhausted budgets return native errors, never successful partial trajectories. No solver runs merely because a block was declared. Each definition reference has schema grid.reactions.reference.v1 or grid.equations.reference.v1 and a name field.

Model quantities are finite numeric values. Reaction initial amounts must also be nonnegative. Unannotated models use the author's consistent numeric coordinates; annotated models use the physical-unit rules below. Currency tags and effectful/volatile calls are not admitted in this profile. Reusable pure numeric functions are described below. The builtin subset is ABS, MIN, MAX, SQRT, EXP, LN, LOG, LOG10, SIN, COS, TAN, ASIN, ACOS, ATAN, SINH, COSH, TANH, POWER (POW), SIGN, FLOOR, CEILING, ROUND, and IF.

Limits are 256 variables, 1,024 directed reaction channels (a reversible declaration counts as two), 1,024 terms per reaction side, 2,048 clauses, 64 expression/dependency levels, and 50,000 expression work units and generated definition nodes. Existing compiler, callable and solver budgets also apply. The owning AST retains every clause and expression; source-free rendering uses assignment spelling, while source formatting preserves the original header and comments. Renaming dynamics-local declarations is refused; state labels and scenario keys need an explicit public rename contract.

Named reusable solutions and observations

REACTION_SOLUTION and EQUATION_SOLUTION take the same definition, time span, and literal scenario/solver options as SOLVE. They request one integration and return a portable grid.dynamics.solution.v1 record. Reading several fields of that bound solution reuses the result through ordinary reactive dependencies; it does not integrate separately for each field. Native solve caching may also reuse an identical request. Declaring or reading the model definition still does no work.

conversion = REACTIONS
  PARAMETER rate = 0.1
  SPECIES substrate INITIAL 10
  SPECIES product INITIAL 0
  LET flux = rate * substrate
 
  REACTION transfer: substrate -> product AT flux
  OBSERVE total = substrate + product
  OBSERVE activity = flux
END
 
run = REACTION_SOLUTION(conversion, 0, 20, {samples: 101})
produced = run.final.product
activity = run.series.activity
work = run.report.work

OBSERVE name = expression exposes a finite numeric derived output. Expressions can read states/species, parameters, TIME, and block-local LET values. Output names share the model namespace; duplicates, shadowing and references to other observations are rejected. Observations are evaluated at returned sample times, not at every solver step, and do not enter the error-control norm. DAE samples use the projected consistent initial state. An invalid/nonfinite observation fails the whole named solution. SOLVE/FINAL/REPORT do not evaluate observations.

The record fields are:

ODE reports use grid.ode.report.v1. DAE reports use grid.dae.report.v1, with named work counters, terminal state, accepted-step range, and initial_consistency (initial/projected scaled residuals, projection iterations, algebraic-variable count and projected state). Report state vectors contain integrated/algebraic variables only. Names preserve authored spelling; record field access uses those names. The outer schema is common to both solver kinds.

At most 256 observations and 2,000,000 numeric trajectory entries are admitted. The solution always returns a record; omit format on these two authoring calls. Both scenario records and all native solver controls remain supported. Solver and output failures never publish a successful partial solution. Backward and zero-length spans retain the native sampling contract.

The native lowering uses ODE_SOLUTION(..., {format: "record"}) or DAE_SOLVE(..., {format: "record"}), then DYNAMICS.SOLUTION to assemble labels and observations through ordinary compiled callbacks. The native record contains trajectory and report. Existing native default numeric results are unchanged. Existing trajectory-query functions can consume run.trajectory; queries interpolate stored samples, including observation columns, rather than invoking the original integrator or re-evaluating observation expressions.

Shared pure numeric LET callbacks now compile their local values once per callback invocation. The compiler retains eager binding order, finite guards, transitive state dependencies, and exact Jacobian/JVP chain rules. Unsupported expressions and unused locals whose errors may be ignored by general LET retain the existing compiled fallback. Dynamics guards deliberately check every declared LET and reaction flux, including unused or zero-net values.

Runnable examples

• Enzyme binding and catalysis: reversible mass action, shared turnover, enzyme/substrate conservation, and an alternative initial enzyme concentration. This retains the full elementary mechanism, without imposing the Michaelis–Menten quasi-steady approximation. Model background: MIT Systems Biology notes.
• Closed population compartments: named incidence and recovery fluxes, conserved population, and changed transmission. Parameters are illustrative, with homogeneous mixing, no births or deaths, and permanent removal from susceptibility. Model background: Duke's SIR model.
• Damped oscillator: adjustable inertia, named spring/total forces, and initial-displacement scenarios. The linear model can be checked against its underdamped analytic solution. Model background: OpenStax damped oscillations.
• RC circuit: parameterized capacitance with an algebraic current constraint and native initial projection. The voltage is supply*(1-EXP(-TIME/(resistance*capacitance))) for the zero initial voltage used here.

These four examples use the author's consistent numeric units. Conservation/analytic checks assess the authored models and numerics separately from successful solver completion. Piecewise forcing via IF is possible, but discontinuities should be split into explicit intervals; these blocks do not promise event-triggered state resets or hybrid scheduling. Imported model definitions/computed definition selection, indexed state families, stochastic kinetics, state-dependent mass, and higher-index DAE reduction remain outside this profile. Existing native trajectory-query functions can analyze sampled results; named threshold events are described below.

Named threshold events

Declare a zero-crossing condition with a name, then request its first crossing:

conversion = REACTIONS
  PARAMETER speed = 0.1
  SPECIES substrate INITIAL 10
  SPECIES product INITIAL 0
  LET flux = speed * substrate
 
  REACTION transfer: substrate -> product AT flux
  EVENT half_converted WHEN substrate = 5 FALLING
  OBSERVE total = substrate + product
END
 
hit = REACTION_EVENT(conversion, "half_converted", 0, 20)
half_life = hit.time
produced = hit.state.product
material = hit.observations.total

EQUATION_EVENT has the same arguments and result. The event condition must be left = right; the solver locates a zero of left - right. RISING, FALLING, or optional ANY filters the residual's direction with respect to increasing time, including backward integration. An exactly zero initial residual counts as a hit. Event names share the state/parameter/observation namespace; there are at most 256 events per definition. Event expressions can use parameters, states, TIME, and shared LET values. They cannot depend on observations or other events.

The result has schema grid.dynamics.event.v1, event, numeric time, named state, and named observations. It is an ordinary reusable value. Observations are evaluated once at the located state; invalid observations fail the whole result. No crossing returns #N/A; invalid event evaluations and exhausted solver budgets remain errors. SOLVE, FINAL, REPORT and SOLUTION do not evaluate event conditions or change their behavior because events were declared.

The optional literal record supports parameter/initial overrides and native integration controls. Omit samples, format, and event_direction: an event returns one hit and uses its declared direction. Reduce max_step when a model may cross more than once in a step. Detection brackets roots at accepted-step endpoints; tangencies and multiple crossings inside a step are not guaranteed. Events do not reset state or continue a hybrid simulation.

For constant-mass index-1 DAEs, initial algebraic projection precedes the first event evaluation. DAE_EVENT uses the existing adaptive implicit integrator and locates a bracketed root with bounded constrained substeps, preserving algebraic equations at the returned state. It does not interpolate algebraic variables between output samples. Its native form is DAE_EVENT(rhs, mass, initial, start, end, event[, options]); both callbacks take (time, state), and its result is one [time, state...] row.

Physical units

Annotate a model's quantities with the existing IS unit:... spelling. A clock declaration specifies the numerical time coordinate:

motion = EQUATIONS
  TIME IS unit:min
  PARAMETER speed IS unit:m/s = 2
  STATE distance IS unit:cm INITIAL 1m
 
  DERIVATIVE(distance) = speed
  EVENT arrived WHEN distance = 5m RISING
  OBSERVE meters IS unit:m = distance
END
 
run = EQUATION_SOLUTION(motion, 0, 1min, {max_step: 5s, atol: 0.1mm})
distance_in_cm = run.final.distance
distance_in_m = run.final.meters

The solver uses the declared coordinates: centimetres and minutes here. Dimensional checking and conversion factors are compiled into ordinary numeric expressions. The example's derivative is therefore 12,000 centimetres per minute. There is no unit parsing inside numerical callbacks.

IS unit:... is available on STATE, SPECIES, ALGEBRAIC, PARAMETER, LET, and OBSERVE; IS number declares a dimensionless quantity. Any annotation or TIME IS selects dimensional admission for the block. Its default clock is seconds. Quantities without annotations infer units from their expressions; bare numbers infer dimensionless coordinates.

Bare numbers at declaration and call boundaries use the target coordinate: INITIAL 10 means ten centimetres for a state declared in centimetres. Explicit quantities must be compatible. Inside governing expressions, units must balance: distance = 5m is valid, while distance = 5 is refused. Bare zero can appear in any dimension, as in an algebraic residual. Arithmetic, comparisons, and conditional branches convert compatible quantities before combining them. Transcendental functions require dimensionless arguments.

Reaction rates remain complete fluxes, with units of amount per model time. MASS_ACTION derives the required dimensions from its reactants. Mixed compatible species coordinates, such as millimoles and moles, receive the appropriate stoichiometric conversion factors. Amount and concentration cannot be mixed implicitly.

Captured inputs and scenario overrides preserve their declared unit metadata. Start/end times, first_step, and max_step accept compatible time quantities. Fixed duration literals such as 1s, 5min, and 2h are supported; calendar months/years are refused. A time-valued reactive input must be numeric and tagged, for example input elapsed IS unit:s = 10; a duration-valued string cell is not a numeric solver coordinate.

An absolute tolerance may be a quantity compatible with every state or an array with one compatible quantity per state. A bare scalar atol uses that same number in each state's declared coordinate. Relative tolerance, sample count, step count and Jacobian bandwidth are dimensionless.

Named solutions and event results retain numeric values in those coordinates and add units.time, units.state.<name>, and units.observations.<name>. Each descriptor contains symbol in normalized Grid type-tag spelling, canonical dimension, and base-unit scale. These are nested records, so both Grid field access and exported JSON retain their named fields. Existing SOLVE, FINAL and REPORT result shapes remain unchanged.

Physical captures unwrap numeric type tags before coordinate conversion and integration. Plain text, booleans, blanks, nonnumeric tagged payloads, and nonfinite values fail the numeric capture check; a tag or unit conversion does not coerce its payload into a number. External model values are captured once per solve; scenario, time, and numeric solver-option inputs are checked at their operation boundary.

Model coordinates must have zero offset: use Kelvin for a temperature state. Celsius/Fahrenheit literals and explicitly tagged inputs are converted once at capture. Custom UNIT declarations use Grid's existing registry.

The physical dynamics example combines chemistry in millimoles/minutes, dimensionless population growth with a physical rate, and an electrical DAE in millivolts/milliseconds.

Reusable numeric functions

Use ordinary Grid function definitions for rate laws, forces, observations and event expressions:

saturating_flux(amount, maximum, half_saturation) = maximum * amount / (half_saturation + amount)
 
conversion = REACTIONS
  PARAMETER maximum = 0.2
  PARAMETER half_saturation = 0.5
  SPECIES substrate INITIAL 1
  SPECIES product INITIAL 0
 
  REACTION transfer: substrate -> product AT saturating_flux(substrate, maximum, half_saturation)
  OBSERVE total = substrate + product
END
 
run = REACTION_SOLUTION(conversion, 0, 10)

Declare helpers before the model. Helpers can call other admitted helpers and use numeric LET locals. Each argument is evaluated once using the caller's scope. Free names in a helper retain its defining module scope; a model parameter with the same spelling does not replace that capture. Captured inputs remain reactive and receive the model's finite numeric checks once per selected operation, including captures used only in a conditional branch. Imported helpers retain their defining callable identities, including private helper calls, through the existing module resolver.

Admission specializes the existing function body to ordinary numeric LET expressions. Units are inferred from each call's argument quantities and checked through the helper body. No second evaluator or function language is introduced. The existing compiled fallback preserves general LET behavior, including ignored errors in unused locals; finite checks still apply to model quantities and returned fluxes or observations.

Effects, unresolved references, recursion, arrays and higher-order callbacks are refused. A helper takes at most 64 arguments, calls nest at most 32 levels, and the existing expression and generated-node budgets apply. A helper used as a derivative coefficient must depend only on solve-constant parameters.

The reusable dynamics examples cover saturating kinetics, a population rate with a lexical capture, and shared mechanical forces.

Coupled model components

COMPONENT includes a preceding model's states and governing equations in the parent's single solve. Qualified names connect components:

vessel = REACTIONS
  PARAMETER leakage = 0.1
  SPECIES material INITIAL 1
  REACTION loss: material -> 0 AT leakage * material
END
 
connected = REACTIONS
  PARAMETER exchange = 0.5
  COMPONENT donor = vessel WITH {parameters: {leakage: 0}}
  COMPONENT receiver = vessel WITH {parameters: {leakage: 0}, initial: {material: 0}}
  REACTION transfer: donor.material -> receiver.material AT exchange * donor.material
  OBSERVE total = donor.material + receiver.material
END
 
run = REACTION_SOLUTION(connected, 0, 2)
remaining = run.final["donor.material"]

Equation components expose changing forces, inflows, or signals with INPUT. Unlike a PARAMETER, an input formula is evaluated at the current model time and state on each callback. Give it a default, then replace that default with an explicit component binding:

body = EQUATIONS
  PARAMETER inertia = 1
  INPUT force = 0
  STATE position INITIAL 0
  STATE velocity INITIAL 0
  DERIVATIVE(position) = velocity
  inertia * DERIVATIVE(velocity) = force
END
 
coupled = EQUATIONS
  COMPONENT first = body WITH {initial: {position: 1}, inputs: {force: second.position - first.position}}
  COMPONENT second = body WITH {inputs: {force: first.position - second.position}}
END
 
motion = EQUATION_SOLUTION(coupled, 0, 1)

Place INPUT and LET before governing clauses. Inputs can use parameters, states, time, and earlier numeric locals. They cannot serve as constant mass coefficients, even if their default happens to be a constant. Component input bindings can use parent states and numeric locals; state feedback is a coupled differential system, not a cycle of computed input definitions.

Components retain declaration order, lexical captures, inherited observations, and named events. Their public labels are qualified strings such as "first.position" or "donor.half". Nested scenario keys use those same strings: {parameters: {"nested.exchange": 2}}. Reading a definition remains passive. An explicit operation integrates the entire composed state vector once.

For physical models, annotate both parent and components. A parent without an explicit TIME uses seconds when physical components select dimensional admission. Child states retain their declared coordinates, child clocks are converted to parent time, and input bindings must have compatible dimensions. Constant mass matrices combine into the same existing ODE/DAE problem.

Composition is bounded by the aggregate 256-state, 256-observation, 256-event, 1,024-reaction-channel and expression limits, with at most 16 component levels. Components must use the parent's notation kind. Model selection uses preceding module definitions or immutable aliases; imported model definitions and computed model selection remain unsupported.

The composed dynamics example combines connected vessels, coupled mechanical bodies, and electrical DAEs with different clock units.

Indexed states, species, and reactions

Use fixed domains to describe repeated compartments or populations. Admission expands each family into ordinary scalar solver coordinates:

compartments = [1, 2, 3]
links = [1, 2]
 
chain = REACTIONS
  PARAMETER transfer_rate = 0.5
  SPECIES amount[i IN compartments] INITIAL 1
 
  REACTION transfer[i IN links]:
    amount[i] -> amount[i + 1] AT transfer_rate * amount[i]
  OBSERVE total = SUM(i IN compartments, amount[i])
END
 
run = REACTION_SOLUTION(chain, 0, 2)
first_amount = run.final["amount[1]"]

STATE, ALGEBRAIC, and named equation rows use the same binders:

regions = ["north", "south"]
 
colonies = EQUATIONS
  PARAMETER growth_rate = 0.1
  STATE population[region IN regions, cohort IN [1, 2]] INITIAL 10 * cohort
 
  EQUATION growth[region IN regions, cohort IN [1, 2]]:
    DERIVATIVE(population[region, cohort]) = growth_rate * population[region, cohort]
  OBSERVE total = SUM(region IN regions, SUM(cohort IN [1, 2], population[region, cohort]))
END

Domains are literal row/column arrays or preceding immutable constants. Each axis contains unique integer, text, or symbol keys of one representation. Text and symbol keys are case-sensitive. Binders are lexical; integer index arithmetic must produce an exact key in the target family. Bare family reads, missing keys, mismatched subscripts, input-driven shapes, and later domain redefinitions are refused. This bounded SUM spelling operates inside indexed dynamics expressions.

Declaration order and domain order determine output columns; the last axis varies fastest. Families can have up to eight axes, with at most 1,024 keys per axis, while the whole model still permits at most 256 expanded states and 1,024 directed reaction channels. The Cartesian product is checked before expansion. Components qualify generated labels, and scenario records use those exact strings, such as "vessel.amount[1]" or "population[\"north\",2]".

The indexed dynamics examples cover a reaction chain, population cohorts, and parallel electrical DAE branches.

Model invariants

An invariant states a condition that every admitted numerical model state must satisfy within an explicit absolute tolerance:

conversion = REACTIONS
  PARAMETER rate = 0.5
  SPECIES substrate INITIAL 1
  SPECIES product INITIAL 0
 
  REACTION transfer: substrate -> product AT rate * substrate
  INVARIANT material: substrate + product = 1 TOLERANCE 1e-8
  INVARIANT nonnegative: substrate >= 0 TOLERANCE 1e-10
END
 
run = REACTION_SOLUTION(conversion, 0, 10)

Relations support =, <=, and >=. Equality checks the absolute difference; inequalities allow at most the stated excess in the forbidden direction. Tolerance must be finite, nonnegative, and constant during the solve. Use a literal or captured PARAMETER; state, time, and dynamic INPUT values cannot set a changing tolerance. A physical tolerance must match the relation's dimension. Bare numbers use the relation's numeric coordinate.

Invariants can use indexed binders, shared numeric helpers, and model locals. Components inherit checks under qualified names. At most 256 expanded checks are admitted across a composed model.

Every solve, final, report, solution, and event operation applies the checks. They run on the initial state after any DAE projection, accepted integration states, requested output samples, and a returned event state. They do not run on rejected trial stages or Newton guesses. Event operations stop at the first requested crossing and do not check unused states beyond it. Reading a model definition still performs no integration or invariant evaluation.

A failed or nonfinite check fails the entire operation with its authored name and time. No successful partial solution is returned. Invariants neither clamp states nor adjust tolerances, and sampled numerical checks do not establish a continuous-time proof between accepted states. Ordinary models without invariants retain the solver's static no-check specialization.

The invariant examples add material conservation, population nonnegativity, and Kirchhoff constraints to the indexed chemistry, biology, and engineering models.