Count-conserving bimolecular CRNs and their mass-action jump chain (CRN-1)

A chemical reaction network (CRN) on a finite set S of species is a finite set of reactions. We consider the count-conserving bimolecular ones, A + B → C + D: two reactant molecules are replaced by two product molecules. The total number n of molecules is then conserved, and the state space is the finite set Counts S n of count vectors.

Under stochastic mass-action kinetics with a common rate constant k, reaction r fires at rate aᵣ(x) = k · ∏ₐ C(xₐ, νₐ), where ν counts the reactants of each species: k·#A·#B for A + B with A ≠ B and k·C(#A, 2) for A + A, that is k times the number of unordered pairs of molecules that can react (Gillespie's combinatorial form, as in Doty, SODA 2014, §2). The continuous-time Markov chain on counts (Anderson–Kurtz 2011) leaves x at rate a₀(x) = ∑ᵣ aᵣ(x); its jump chain fires r with probability aᵣ(x) / a₀(x). A terminal state (a₀(x) = 0) is made absorbing.

structure Crn.Reaction (S : Type u_1) : Type u_1

A count-conserving bimolecular reaction A + B → C + D (roadmap CRN-1): the unordered pair of reactant species s(A, B) is replaced by the unordered pair of product species s(C, D). Equal species (A = B or C = D) are allowed.

• reactants : Sym2 S

  The two reactant species.

• products : Sym2 S

  The two product species.

def Crn.instDecidableEqReaction.decEq {S✝ : Type u_1} [DecidableEq S✝] (x✝ x✝¹ : Reaction S✝) : Decidable (x✝ = x✝¹)

Equations

• Crn.instDecidableEqReaction.decEq { reactants := a, products := a_1 } { reactants := b, products := b_1 } = if h : a = b then h ▸ if h : a_1 = b_1 then h ▸ isTrue ⋯ else isFalse ⋯ else isFalse ⋯

@[implicit_reducible]

instance Crn.instDecidableEqReaction {S✝ : Type u_1} [DecidableEq S✝] : DecidableEq (Reaction S✝)

Equations

• Crn.instDecidableEqReaction = Crn.instDecidableEqReaction.decEq

structure Crn.Network (S : Type u_1) : Type u_1

A count-conserving bimolecular CRN (roadmap CRN-1): a nonempty finite set of reactions, none of which leaves the counts unchanged. All reactions share one rate constant, which is a parameter of the kinetics (Network.jumpKernel).

• reactions : Finset (Reaction S)

  The reactions.

• nonempty : self.reactions.Nonempty

  There is at least one reaction.

• nontrivial (r : Reaction S) : r ∈ self.reactions → r.reactants ≠ r.products

  No reaction is a no-op.

@[reducible, inline]

abbrev Crn.Counts (S : Type u_2) [Fintype S] (n : ℕ) : Type u_2

Count vectors of n molecules: the number of molecules of each species, summing to n.

Equations

• Crn.Counts S n = { x : S → ℕ // ∑ a : S, x a = n }

@[implicit_reducible]

instance Crn.Counts.instFintype {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} : Fintype (Counts S n)

Count vectors of n molecules form a finite type, enumerated by Finset.piAntidiag.

Equations

• Crn.Counts.instFintype = Fintype.subtype (Finset.univ.piAntidiag n) ⋯

def Crn.Reaction.consumed {S : Type u_1} [DecidableEq S] (r : Reaction S) (a : S) : ℕ

Number of reactant molecules of species a (0, 1 or 2).

Equations

• r.consumed a = Multiset.count a r.reactants.toMultiset

def Crn.Reaction.produced {S : Type u_1} [DecidableEq S] (r : Reaction S) (a : S) : ℕ

Number of product molecules of species a (0, 1 or 2).

Equations

• r.produced a = Multiset.count a r.products.toMultiset

theorem Crn.Reaction.sum_consumed {S : Type u_1} [DecidableEq S] [Fintype S] (r : Reaction S) : ∑ a : S, r.consumed a = 2

A reaction consumes two molecules.

theorem Crn.Reaction.sum_produced {S : Type u_1} [DecidableEq S] [Fintype S] (r : Reaction S) : ∑ a : S, r.produced a = 2

A reaction produces two molecules.

noncomputable def Crn.Reaction.propensity {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} (k : ℝ) (r : Reaction S) (x : Counts S n) : ℝ

Stochastic mass-action propensity of r at counts x with rate constant k, in Gillespie's combinatorial form k · ∏ₐ C(xₐ, νₐ), ν the reactant multiplicities (propensity_of_ne, propensity_of_eq). It vanishes iff r lacks reactant molecules.

Equations

• Crn.Reaction.propensity k r x = k * ∏ a : S, ↑((↑x a).choose (r.consumed a))

theorem Crn.Reaction.count_toMultiset_mk {S : Type u_1} [DecidableEq S] (A B a : S) : Multiset.count a s(A, B).toMultiset = (if a = A then 1 else 0) + if a = B then 1 else 0

Multiplicity of the species a in the unordered pair s(A, B).

theorem Crn.Reaction.consumed_le_of_propensity_ne_zero {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} {k : ℝ} {r : Reaction S} {x : Counts S n} (h : propensity k r x ≠ 0) (a : S) : r.consumed a ≤ ↑x a

A reaction with nonzero propensity has all its reactant molecules.

theorem Crn.Reaction.propensity_nonneg {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} {k : ℝ} (hk : 0 ≤ k) (r : Reaction S) (x : Counts S n) : 0 ≤ propensity k r x

Propensities are nonnegative for a nonnegative rate constant.

theorem Crn.Reaction.propensity_of_ne {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} (k : ℝ) (r : Reaction S) (x : Counts S n) {A B : S} (hr : r.reactants = s(A, B)) (hAB : A ≠ B) : propensity k r x = k * ↑(↑x A) * ↑(↑x B)

For A + B → ⋯ with A ≠ B the propensity is k·#A·#B.

theorem Crn.Reaction.propensity_of_eq {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} (k : ℝ) (r : Reaction S) (x : Counts S n) {A : S} (hr : r.reactants = s(A, A)) : propensity k r x = k * ↑((↑x A).choose 2)

For A + A → ⋯ the propensity is k·C(#A, 2).

def Crn.Counts.react {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} (x : Counts S n) (r : Reaction S) : Counts S n

The counts after one firing of r: xₐ - νₐ + ν'ₐ with ν, ν' the reactant and product multiplicities. If r lacks reactant molecules, x is returned unchanged (such a firing has probability zero in every chain below).

Equations

• x.react r = if h : ∀ (a : S), r.consumed a ≤ ↑x a then ⟨fun (a : S) => ↑x a - r.consumed a + r.produced a, ⋯⟩ else x

theorem Crn.Counts.react_val {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} {x : Counts S n} {r : Reaction S} (h : ∀ (a : S), r.consumed a ≤ ↑x a) (a : S) : ↑(x.react r) a = ↑x a - r.consumed a + r.produced a

The counts after firing an applicable reaction.

theorem Crn.Counts.react_ne_self {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} {x : Counts S n} {r : Reaction S} (hr : r.reactants ≠ r.products) (h : ∀ (a : S), r.consumed a ≤ ↑x a) : x.react r ≠ x

Firing an applicable reaction whose products differ from its reactants changes the counts.

noncomputable def Crn.Network.totalPropensity {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} (N : Network S) (k : ℝ) (x : Counts S n) : ℝ

Total propensity a₀(x) = ∑ᵣ aᵣ(x), the rate at which the counts leave x.

Equations

• N.totalPropensity k x = ∑ r ∈ N.reactions, Crn.Reaction.propensity k r x

noncomputable def Crn.Network.nextReaction {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} (N : Network S) {k : ℝ} (hk : 0 < k) (x : Counts S n) (h : N.totalPropensity k x ≠ 0) : Dynamics.Distribution ↥N.reactions

The reaction fired by the jump chain from x: r with probability aᵣ(x) / a₀(x).

Equations

• N.nextReaction hk x h = Crn.normalize (fun (r : ↥N.reactions) => Crn.Reaction.propensity k (↑r) x) ⋯ ⋯

noncomputable def Crn.Network.jumpKernel {S : Type u_1} [Fintype S] [DecidableEq S] (N : Network S) (k : ℝ) (hk : 0 < k) (n : ℕ) : Dynamics.Kernel (Counts S n)

The jump chain of stochastic mass-action kinetics with common rate constant k > 0 (roadmap CRN-1): from x, fire reaction r with probability aᵣ(x) / a₀(x); a terminal state (a₀(x) = 0) is absorbing.

Equations

• N.jumpKernel k hk n x = if h : N.totalPropensity k x = 0 then Dynamics.Distribution.point x else (N.nextReaction hk x h).map fun (r : ↥N.reactions) => x.react ↑r

theorem Crn.Network.totalPropensity_mul_expect {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} (N : Network S) {k : ℝ} (hk : 0 < k) (x : Counts S n) (f : Counts S n → ℝ) : N.totalPropensity k x * (N.jumpKernel k hk n x).expect f = ∑ r ∈ N.reactions, Reaction.propensity k r x * f (x.react r)

The jump chain's expectation times the total propensity, valid also at terminal states: a₀(x)·E[f] = ∑ᵣ aᵣ(x)·f(x.react r).

theorem Crn.Network.jumpKernel_weight_self {S : Type u_1} [Fintype S] [DecidableEq S] {n : ℕ} (N : Network S) {k : ℝ} (hk : 0 < k) (x : Counts S n) (h : N.totalPropensity k x ≠ 0) : (N.jumpKernel k hk n x).weight x = 0

The jump chain really jumps: from a non-terminal state it never stays put (every reaction of the network changes the counts).