Problems with reactions where two monomers form another monomer VS two monomers form a dimer

[Quim98]

Hello, I have found that when defining reactions where two monomers of the same molecule bind to form either a dimer or another monomer, the model does not arrive to the mathematically-derived steady-state. For example, if we define a reaction where two monomers form a dimer, the system only gets to the mathematicall-derived steady state if the dimer concentration is divided by 2. In contrast, if we define a reaction where two monomers form another monomer, the system only gets to the mathematically-derived steady state if the backward reaction rate is divided by two. So, both of the explained reaction examples can be defined by a single ODE, where D is the dimer, M is the free monomer, k_B is the binding rate and k_U is the unbinding rate:

$\frac{dD}{dt}=k_BM^2-k_UD$ (1)

Taking into account that the total monomer (free and bounded) is constant M_0=M+2D, and that k_D=k_U/k_B, the steady-state is found via the next expression:

$M*=\frac{\frac{-k_D}{2}+\sqrt{\frac{k_D}{4}+2k_DM_0}}{2}$ (2)

Given k_B=100, k_U=0.01 and M_0=1e-4, then the free monomer (M*) at the steady state is found at 5e-5

Beggining with the examples, the first model defines a reaction where two monomers of the same molecule bind to form a dimer, where the steady states of the dimer and monomer form are calculated via the equation 2 and used as initial conditions. If the result is printed we can see that the system arrives to the mathematically-derived steady-state for the free monomer if the dimer concentration at the steady state is divided by two, as the package seems to be mutiplying by two the inputed value for the initial conditions.

from pysb import Model, Parameter, Compartment, Monomer, Parameter, Rule, Initial, Observable
from pysb.simulator import ScipyOdeSimulator
import numpy as np
Model()
Monomer('M', ['b'])
Parameter('K_b', 1E+2)
Parameter('K_u', 1E-2)
M_free = M(b=None)
D_free = M(b=1) % M(b=1)
Rule('IL_binding', M_free + M_free | D_free, *[K_b, K_u])
Kd = 1E-2/1E+2
M_eq = (-Kd/2 + (Kd**2/4 + 2*Kd*1e-4)**(1/2))/2
D_eq = M_eq**2/Kd
Parameter('M_0', M_eq)
Initial(M_free, M_0)
Parameter('D_0', D_eq/2) # Initial condition of the dimer divided by 2
Initial(D_free, D_0)
Observable('M_free', M_free)
Observable('D_free', D_free)
t = np.arange(0,100,0.001)
simulator = ScipyOdeSimulator(model, tspan=t, compiler='python').run()
result = simulator.all
print("Calculated monomer at steady state: "+str(M_eq))
print("Calculated dimer at steady state: "+str(D_eq))
print("Free monomer simulation: "+str(result['M_free']))
print("Free dimer simulation: "+str(result['D_free']))

Instead, we can define the dimer species as another monomer. As in the last example, the steady states of the dimer and monomer form are calculated via the equation 2 and used as initial conditions. If the result is printed we can see that the system arrives to the mathematically-derived steady-state for the free monomer if the unbinding rate is divided by two, as the package seems to be mutiplying by two the inputed value for the unbinding rate.

from pysb import Model, Parameter, Compartment, Monomer, Parameter, Rule, Initial, Observable
from pysb.simulator import ScipyOdeSimulator
import numpy as np
Model()
Monomer('M')
Monomer('D')
Parameter('K_b', 1E+2)
Parameter('K_u', 1E-2/2) # Unbinding rate divided by 2
M_free = M()
D_free = D()
Rule('IL_binding', M_free + M_free | D_free, *[K_b, K_u])
Kd = 1E-2/1E+2  # To calculate the steady state the unbinding rate is not divided by 2
M_eq = (-Kd/2 + (Kd**2/4 + 2*Kd*1e-4)**(1/2))/2
D_eq = M_eq**2/Kd
Parameter('M_0', M_eq)
Initial(M_free, M_0)
Parameter('D_0', D_eq)
Initial(D_free, D_0)
Observable('M_free', M_free)
Observable('D_free', D_free)
t = np.arange(0,100,0.001)
simulator = ScipyOdeSimulator(model, tspan=t, compiler='python').run()
result = simulator.all
print("Calculated monomer at steady state: "+str(M_eq))
print("Calculated dimer at steady state: "+str(D_eq))
print("Free monomer simulation: "+str(result['M_free']))
print("Free dimer simulation: "+str(result['D_free']))

These are the only two instances that I've found this strange behaviour, but I have no clue where this problem could source from, as I have reviewed multiple times the mathematically-derived steady-state and it seems to be correct.

Thanks in advance