Odd behavior when modeling second order reaction with one species

[jadeshi]

Hi,
I'm running into a very strange issue when trying to model a second order reaction that contains only one species, wondering if I could get some help. Specifically, I'm comparing the following reactions:

1. A + A -> None
2. A + A -> A

What's odd to me is when I run both reactions and plot concentration over time I get exactly overlapping curves. I used the following code:

from pysb import *
from pysb.simulator import ScipyOdeSimulator
import matplotlib.pyplot as plt
import numpy as np

def deg_create(tspan):
    Model()
    Monomer('A')
    Monomer('B')
    Parameter('k', 1.3)
    Rule('rxn1', A() + A() >> A(), k)
    Rule('rxn2', B() + B() >> None, k)
    Parameter('A0', 1.1)
    Initial(A(), A0)
    Initial(B(), A0)
    Observable('obsA', A())
    Observable('obsB', B())
    res = ScipyOdeSimulator(model, tspan=tspan, compiler='python').run().all
    return res

dt = 0.001
t = np.arange(0,10,dt)
resa = deg_create(t)['obsA'] # 2 to 0
resb = deg_create(t)['obsB'] # 2 to 1

plt.plot(resb, color='red', label=' 2 to 0',linewidth=5)
plt.plot(resa, color='blue', label='2 to 1')
plt.legend()
plt.xlabel('Time', fontsize=14)
plt.ylabel('Concentration', fontsize=14)
plt.show()

I get the following plot:

This obviously seems incorrect to me, since the two reactions are clearly not the same. The rate of change of the first reaction should in theory be half of the second since you're generating A half as fast as you're depleting it based on the stoichiometry.

Curiously though, I did a follow-up sanity check and added in an intermediate species to turn the second reaction into a two-step reaction:

A + A -> A', k=1.3
A' -> A, k=10000

where in this case the A' is being converted infinitely fast back into A, so it's functionally equivalent to A + A -> A, I get different, seemingly correct kinetics.

I'm wondering, am I defining the reaction incorrectly to begin with? Is it not possible to write a reaction like this with one species? Thanks very much for the help!

[lh64]

Hi there. This actually has to do with how PySB (or more accurately, BioNetGen) handles symmetry in rules. Basically, for every rule, two symmetry factors are calculated: a "multiplicity factor", m, and a "path factor", p. These factors are included in the rate expression of the rule by automatically multiplying by the rate constant k, i.e., the "effective" rate constant is p*m*k.

The multiplicity factor has to do with the combinatorics of the reactant species. For example, in a rule like B() + B() >> None, the multiplicity factor m=1/2. This is because there are X_B*(X_B-1)/2 different ways in which two B molecules can interact, where X_B represents the number of B molecules. For ODE simulations, like yours, X_B*(X_B-1) is approximately X_B^2, so the reaction rate for the rule is calculated as rate = 0.5*k*X_B^2, which gives the familiar reaction rate equation. It is important to recognize that the factor of 0.5 gets automatically included in the rate expression.

The path factor has to do with how many alternate paths exist for converting the reactants in a rule into the products. For example, in a rule like A() + A() >> A(), one of the A() molecules is degraded. Since there are two A() molecules on the reactant side of the rule, that means there are two different paths for converting the reactants into products, i.e., the first molecule could degrade or the second molecule could degrade. Therefore, p=2. Since this rule also has the same multiplicity factor of 1/2, the rate for this rule is rate = 2*0.5*k*X_A^2 = k*X_A^2, i.e., the factors of 2 and 0.5 cancel out.

Finally, the reason the two rules A() + A() >> A() and B() + B() >> None (with identical rate constants k) are producing identical outputs is because of the difference in stoichiometries of the reactant species. In the first rule, A() has a stoichiometry of -1, while in the second rule B() has a stoichiometry of -2. When you put it all together, you get the differential equations dX_A/dt = -k*X_A^2 and dX_B/dt = -2*0.5*k*X_B^2 = -k*X_B^2. So, unintentionally, you've created two reactions with identical rate expressions. This is a consequence of the graph-theoretic nature of the underlying rule-based formalism, which is admittedly sometimes confusing.

BTW, you can see the rate expressions for all reactions, including the multiplicity factors, by printing the reactions to the screen. After the simulation completes, just include in your code:

for r in model.reactions:
print(r)

Each reaction is a dictionary. There should be an entry called rate where you can see the exact expression that's used.

I hope this helps. If have any additional questions, just let us know.

--Leonard

[jadeshi]

Ah, thanks very much for clarifying, I was unaware of the multiplicity and path factors. So hypothetically, if I were to build a model using experimentally determined rate constants, is it possible that the simulation could underestimate the rate by applying an extra multiplicity factor, since experimentally, you're typically limited to measuring the "effective" rate constant with the multiplicity information implicitly accounted for?

[lh64]

Yes, that’s correct. The rate constants in a PySB model are what we call “instance rates”, i.e., the rate at which one instance of a rule occurs. In this case, the rate at which one A() molecule degrades or the rate at which any particular pair of A() or B() molecules react. Experimentally, it’s not always possible to measure these rates, as you said. So, it’s definitely something to be aware of.