Catalyst/DifferentialEquations: Gillespie simulation not working for large chemical reaction networks

[uagaug]

Transferring from this thread on Discourse.

I have a large chemical reaction network that involves protein synthesis, degradation, binding, and catalysis. I was able to construct the network using Catalyst, and deterministic ODE simulation in Julia returns the same result as that in Python. However when I tried to run Gillespie simulation it finished within 1 second, returned all specie counts being 0, and gave no error message.

I tried to increase tspan to 10^12 but the simulation did not finish after 12 hours and the ActivityMonitor on my Mac showed that Julia was using 1.2 TB of memory (the machine only has 128 GB memory). The same reaction network finishes within 3 hours with a hand written Gillespie code in Python.

I attach the code below, just wondering if there are some settings that I might have missed. Thanks!

using DifferentialEquations, Catalyst

rn = @reaction_network begin
	(kon1, koff1), x1 + n11D ↔ x1Vn11D
	(kon1, koff1), x1 + n12D ↔ x1Vn12D
	(kon1, koff1), x2 + n21D ↔ x2Vn21D
	(kon1, koff1), x2 + n22D ↔ x2Vn22D
	(kon2, koff2), y1 + x1Vn11D ↔ y1Vx1Vn11D
	(kon2, koff2), y1 + x1Vn12D ↔ y1Vx1Vn12D
	(kon2, koff2), y1 + x2Vn21D ↔ y1Vx2Vn21D
	(kon2, koff2), y1 + x2Vn22D ↔ y1Vx2Vn22D
	(kon2, koff2), y2 + x1Vn11D ↔ y2Vx1Vn11D
	(kon2, koff2), y2 + x1Vn12D ↔ y2Vx1Vn12D
	(kon2, koff2), y2 + x2Vn21D ↔ y2Vx2Vn21D
	(kon2, koff2), y2 + x2Vn22D ↔ y2Vx2Vn22D
	(kon_p1, koff_p1), y1Vx1Vn11D + n11D ↔ y1Vx1Vn11D_n11D
	(kon_p1, koff_p1), y1Vx1Vn11D + n21D ↔ y1Vx1Vn11D_n21D
	(kon_p1, koff_p1), y1Vx2Vn21D + n11D ↔ y1Vx2Vn21D_n11D
	(kon_p1, koff_p1), y1Vx2Vn21D + n21D ↔ y1Vx2Vn21D_n21D
	(kon_p2, koff_p2), y2Vx1Vn12D + n12D ↔ y2Vx1Vn12D_n12D
	(kon_p2, koff_p2), y2Vx1Vn12D + n22D ↔ y2Vx1Vn12D_n22D
	(kon_p2, koff_p2), y2Vx2Vn22D + n12D ↔ y2Vx2Vn22D_n12D
	(kon_p2, koff_p2), y2Vx2Vn22D + n22D ↔ y2Vx2Vn22D_n22D
	kcat1, y1Vx1Vn11D_n11D → y1Vx1Vn11D + n11
	kcat1, y1Vx1Vn11D_n21D → y1Vx1Vn11D + n21
	kcat1, y1Vx2Vn21D_n11D → y1Vx2Vn21D + n11
	kcat1, y1Vx2Vn21D_n21D → y1Vx2Vn21D + n21
	kcat2, y2Vx1Vn12D_n12D → y2Vx1Vn12D + n12
	kcat2, y2Vx1Vn12D_n22D → y2Vx1Vn12D + n22
	kcat2, y2Vx2Vn22D_n12D → y2Vx2Vn22D + n12
	kcat2, y2Vx2Vn22D_n22D → y2Vx2Vn22D + n22
	(kon_p1, koff_p1), y1Vx1Vn11D + x1Vn11D ↔ y1Vx1Vn11D_x1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11D + x2Vn21D ↔ y1Vx1Vn11D_x2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21D + x1Vn11D ↔ y1Vx2Vn21D_x1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21D + x2Vn21D ↔ y1Vx2Vn21D_x2Vn21D
	(kon_p2, koff_p2), y2Vx1Vn12D + x1Vn12D ↔ y2Vx1Vn12D_x1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12D + x2Vn22D ↔ y2Vx1Vn12D_x2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22D + x1Vn12D ↔ y2Vx2Vn22D_x1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22D + x2Vn22D ↔ y2Vx2Vn22D_x2Vn22D
	kcat1, y1Vx1Vn11D_x1Vn11D → y1Vx1Vn11D + x1Vn11
	kcat1, y1Vx1Vn11D_x2Vn21D → y1Vx1Vn11D + x2Vn21
	kcat1, y1Vx2Vn21D_x1Vn11D → y1Vx2Vn21D + x1Vn11
	kcat1, y1Vx2Vn21D_x2Vn21D → y1Vx2Vn21D + x2Vn21
	kcat2, y2Vx1Vn12D_x1Vn12D → y2Vx1Vn12D + x1Vn12
	kcat2, y2Vx1Vn12D_x2Vn22D → y2Vx1Vn12D + x2Vn22
	kcat2, y2Vx2Vn22D_x1Vn12D → y2Vx2Vn22D + x1Vn12
	kcat2, y2Vx2Vn22D_x2Vn22D → y2Vx2Vn22D + x2Vn22
	(kon_p1, koff_p1), y1Vx1Vn11D + y1Vx1Vn11D ↔ y1Vx1Vn11D_y1Vx1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11D + y1Vx2Vn21D ↔ y1Vx1Vn11D_y1Vx2Vn21D
	(kon_p1, koff_p1), y1Vx1Vn11D + y2Vx1Vn11D ↔ y1Vx1Vn11D_y2Vx1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11D + y2Vx2Vn21D ↔ y1Vx1Vn11D_y2Vx2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21D + y1Vx1Vn11D ↔ y1Vx2Vn21D_y1Vx1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21D + y1Vx2Vn21D ↔ y1Vx2Vn21D_y1Vx2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21D + y2Vx1Vn11D ↔ y1Vx2Vn21D_y2Vx1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21D + y2Vx2Vn21D ↔ y1Vx2Vn21D_y2Vx2Vn21D
	(kon_p2, koff_p2), y2Vx1Vn12D + y1Vx1Vn12D ↔ y2Vx1Vn12D_y1Vx1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12D + y1Vx2Vn22D ↔ y2Vx1Vn12D_y1Vx2Vn22D
	(kon_p2, koff_p2), y2Vx1Vn12D + y2Vx1Vn12D ↔ y2Vx1Vn12D_y2Vx1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12D + y2Vx2Vn22D ↔ y2Vx1Vn12D_y2Vx2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22D + y1Vx1Vn12D ↔ y2Vx2Vn22D_y1Vx1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22D + y1Vx2Vn22D ↔ y2Vx2Vn22D_y1Vx2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22D + y2Vx1Vn12D ↔ y2Vx2Vn22D_y2Vx1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22D + y2Vx2Vn22D ↔ y2Vx2Vn22D_y2Vx2Vn22D
	kcat1, y1Vx1Vn11D_y1Vx1Vn11D → y1Vx1Vn11D + y1Vx1Vn11
	kcat1, y1Vx1Vn11D_y1Vx2Vn21D → y1Vx1Vn11D + y1Vx2Vn21
	kcat1, y1Vx1Vn11D_y2Vx1Vn11D → y1Vx1Vn11D + y2Vx1Vn11
	kcat1, y1Vx1Vn11D_y2Vx2Vn21D → y1Vx1Vn11D + y2Vx2Vn21
	kcat1, y1Vx2Vn21D_y1Vx1Vn11D → y1Vx2Vn21D + y1Vx1Vn11
	kcat1, y1Vx2Vn21D_y1Vx2Vn21D → y1Vx2Vn21D + y1Vx2Vn21
	kcat1, y1Vx2Vn21D_y2Vx1Vn11D → y1Vx2Vn21D + y2Vx1Vn11
	kcat1, y1Vx2Vn21D_y2Vx2Vn21D → y1Vx2Vn21D + y2Vx2Vn21
	kcat2, y2Vx1Vn12D_y1Vx1Vn12D → y2Vx1Vn12D + y1Vx1Vn12
	kcat2, y2Vx1Vn12D_y1Vx2Vn22D → y2Vx1Vn12D + y1Vx2Vn22
	kcat2, y2Vx1Vn12D_y2Vx1Vn12D → y2Vx1Vn12D + y2Vx1Vn12
	kcat2, y2Vx1Vn12D_y2Vx2Vn22D → y2Vx1Vn12D + y2Vx2Vn22
	kcat2, y2Vx2Vn22D_y1Vx1Vn12D → y2Vx2Vn22D + y1Vx1Vn12
	kcat2, y2Vx2Vn22D_y1Vx2Vn22D → y2Vx2Vn22D + y1Vx2Vn22
	kcat2, y2Vx2Vn22D_y2Vx1Vn12D → y2Vx2Vn22D + y2Vx1Vn12
	kcat2, y2Vx2Vn22D_y2Vx2Vn22D → y2Vx2Vn22D + y2Vx2Vn22
	(kon_p1, koff_p1), y1Vx1Vn11D + y2 ↔ y1Vx1Vn11D_y2
	(kon_p1, koff_p1), y1Vx2Vn21D + y2 ↔ y1Vx2Vn21D_y2
	(kon_p2, koff_p2), y2Vx1Vn12D + y1 ↔ y2Vx1Vn12D_y1
	(kon_p2, koff_p2), y2Vx2Vn22D + y1 ↔ y2Vx2Vn22D_y1
	kcat1, y1Vx1Vn11D_y2 → y1Vx1Vn11D + y2D
	kcat1, y1Vx2Vn21D_y2 → y1Vx2Vn21D + y2D
	kcat2, y2Vx1Vn12D_y1 → y2Vx1Vn12D + y1D
	kcat2, y2Vx2Vn22D_y1 → y2Vx2Vn22D + y1D
	(kon1, koff1), x1 + n11 ↔ x1Vn11
	(kon1, koff1), x1 + n12 ↔ x1Vn12
	(kon1, koff1), x2 + n21 ↔ x2Vn21
	(kon1, koff1), x2 + n22 ↔ x2Vn22
	(kon2, koff2), y1 + x1Vn11 ↔ y1Vx1Vn11
	(kon2, koff2), y1 + x1Vn12 ↔ y1Vx1Vn12
	(kon2, koff2), y1 + x2Vn21 ↔ y1Vx2Vn21
	(kon2, koff2), y1 + x2Vn22 ↔ y1Vx2Vn22
	(kon2, koff2), y2 + x1Vn11 ↔ y2Vx1Vn11
	(kon2, koff2), y2 + x1Vn12 ↔ y2Vx1Vn12
	(kon2, koff2), y2 + x2Vn21 ↔ y2Vx2Vn21
	(kon2, koff2), y2 + x2Vn22 ↔ y2Vx2Vn22
	(kon_p1, koff_p1), y1Vx1Vn11 + n11D ↔ y1Vx1Vn11_n11D
	(kon_p1, koff_p1), y1Vx1Vn11 + n21D ↔ y1Vx1Vn11_n21D
	(kon_p1, koff_p1), y1Vx2Vn21 + n11D ↔ y1Vx2Vn21_n11D
	(kon_p1, koff_p1), y1Vx2Vn21 + n21D ↔ y1Vx2Vn21_n21D
	(kon_p2, koff_p2), y2Vx1Vn12 + n12D ↔ y2Vx1Vn12_n12D
	(kon_p2, koff_p2), y2Vx1Vn12 + n22D ↔ y2Vx1Vn12_n22D
	(kon_p2, koff_p2), y2Vx2Vn22 + n12D ↔ y2Vx2Vn22_n12D
	(kon_p2, koff_p2), y2Vx2Vn22 + n22D ↔ y2Vx2Vn22_n22D
	kcat1, y1Vx1Vn11_n11D → y1Vx1Vn11 + n11
	kcat1, y1Vx1Vn11_n21D → y1Vx1Vn11 + n21
	kcat1, y1Vx2Vn21_n11D → y1Vx2Vn21 + n11
	kcat1, y1Vx2Vn21_n21D → y1Vx2Vn21 + n21
	kcat2, y2Vx1Vn12_n12D → y2Vx1Vn12 + n12
	kcat2, y2Vx1Vn12_n22D → y2Vx1Vn12 + n22
	kcat2, y2Vx2Vn22_n12D → y2Vx2Vn22 + n12
	kcat2, y2Vx2Vn22_n22D → y2Vx2Vn22 + n22
	(kon_p1, koff_p1), y1Vx1Vn11 + x1Vn11D ↔ y1Vx1Vn11_x1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11 + x2Vn21D ↔ y1Vx1Vn11_x2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21 + x1Vn11D ↔ y1Vx2Vn21_x1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21 + x2Vn21D ↔ y1Vx2Vn21_x2Vn21D
	(kon_p2, koff_p2), y2Vx1Vn12 + x1Vn12D ↔ y2Vx1Vn12_x1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12 + x2Vn22D ↔ y2Vx1Vn12_x2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22 + x1Vn12D ↔ y2Vx2Vn22_x1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22 + x2Vn22D ↔ y2Vx2Vn22_x2Vn22D
	kcat1, y1Vx1Vn11_x1Vn11D → y1Vx1Vn11 + x1Vn11
	kcat1, y1Vx1Vn11_x2Vn21D → y1Vx1Vn11 + x2Vn21
	kcat1, y1Vx2Vn21_x1Vn11D → y1Vx2Vn21 + x1Vn11
	kcat1, y1Vx2Vn21_x2Vn21D → y1Vx2Vn21 + x2Vn21
	kcat2, y2Vx1Vn12_x1Vn12D → y2Vx1Vn12 + x1Vn12
	kcat2, y2Vx1Vn12_x2Vn22D → y2Vx1Vn12 + x2Vn22
	kcat2, y2Vx2Vn22_x1Vn12D → y2Vx2Vn22 + x1Vn12
	kcat2, y2Vx2Vn22_x2Vn22D → y2Vx2Vn22 + x2Vn22
	(kon_p1, koff_p1), y1Vx1Vn11 + y1Vx1Vn11D ↔ y1Vx1Vn11_y1Vx1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11 + y1Vx2Vn21D ↔ y1Vx1Vn11_y1Vx2Vn21D
	(kon_p1, koff_p1), y1Vx1Vn11 + y2Vx1Vn11D ↔ y1Vx1Vn11_y2Vx1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11 + y2Vx2Vn21D ↔ y1Vx1Vn11_y2Vx2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21 + y1Vx1Vn11D ↔ y1Vx2Vn21_y1Vx1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21 + y1Vx2Vn21D ↔ y1Vx2Vn21_y1Vx2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21 + y2Vx1Vn11D ↔ y1Vx2Vn21_y2Vx1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21 + y2Vx2Vn21D ↔ y1Vx2Vn21_y2Vx2Vn21D
	(kon_p2, koff_p2), y2Vx1Vn12 + y1Vx1Vn12D ↔ y2Vx1Vn12_y1Vx1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12 + y1Vx2Vn22D ↔ y2Vx1Vn12_y1Vx2Vn22D
	(kon_p2, koff_p2), y2Vx1Vn12 + y2Vx1Vn12D ↔ y2Vx1Vn12_y2Vx1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12 + y2Vx2Vn22D ↔ y2Vx1Vn12_y2Vx2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22 + y1Vx1Vn12D ↔ y2Vx2Vn22_y1Vx1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22 + y1Vx2Vn22D ↔ y2Vx2Vn22_y1Vx2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22 + y2Vx1Vn12D ↔ y2Vx2Vn22_y2Vx1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22 + y2Vx2Vn22D ↔ y2Vx2Vn22_y2Vx2Vn22D
	kcat1, y1Vx1Vn11_y1Vx1Vn11D → y1Vx1Vn11 + y1Vx1Vn11
	kcat1, y1Vx1Vn11_y1Vx2Vn21D → y1Vx1Vn11 + y1Vx2Vn21
	kcat1, y1Vx1Vn11_y2Vx1Vn11D → y1Vx1Vn11 + y2Vx1Vn11
	kcat1, y1Vx1Vn11_y2Vx2Vn21D → y1Vx1Vn11 + y2Vx2Vn21
	kcat1, y1Vx2Vn21_y1Vx1Vn11D → y1Vx2Vn21 + y1Vx1Vn11
	kcat1, y1Vx2Vn21_y1Vx2Vn21D → y1Vx2Vn21 + y1Vx2Vn21
	kcat1, y1Vx2Vn21_y2Vx1Vn11D → y1Vx2Vn21 + y2Vx1Vn11
	kcat1, y1Vx2Vn21_y2Vx2Vn21D → y1Vx2Vn21 + y2Vx2Vn21
	kcat2, y2Vx1Vn12_y1Vx1Vn12D → y2Vx1Vn12 + y1Vx1Vn12
	kcat2, y2Vx1Vn12_y1Vx2Vn22D → y2Vx1Vn12 + y1Vx2Vn22
	kcat2, y2Vx1Vn12_y2Vx1Vn12D → y2Vx1Vn12 + y2Vx1Vn12
	kcat2, y2Vx1Vn12_y2Vx2Vn22D → y2Vx1Vn12 + y2Vx2Vn22
	kcat2, y2Vx2Vn22_y1Vx1Vn12D → y2Vx2Vn22 + y1Vx1Vn12
	kcat2, y2Vx2Vn22_y1Vx2Vn22D → y2Vx2Vn22 + y1Vx2Vn22
	kcat2, y2Vx2Vn22_y2Vx1Vn12D → y2Vx2Vn22 + y2Vx1Vn12
	kcat2, y2Vx2Vn22_y2Vx2Vn22D → y2Vx2Vn22 + y2Vx2Vn22
	(kon_p1, koff_p1), y1Vx1Vn11 + y2 ↔ y1Vx1Vn11_y2
	(kon_p1, koff_p1), y1Vx2Vn21 + y2 ↔ y1Vx2Vn21_y2
	(kon_p2, koff_p2), y2Vx1Vn12 + y1 ↔ y2Vx1Vn12_y1
	(kon_p2, koff_p2), y2Vx2Vn22 + y1 ↔ y2Vx2Vn22_y1
	kcat1, y1Vx1Vn11_y2 → y1Vx1Vn11 + y2D
	kcat1, y1Vx2Vn21_y2 → y1Vx2Vn21 + y2D
	kcat2, y2Vx1Vn12_y1 → y2Vx1Vn12 + y1D
	kcat2, y2Vx2Vn22_y1 → y2Vx2Vn22 + y1D
	(kon2, koff2), y2D + x1Vn11D ↔ y2DVx1Vn11D
	(kon2, koff2), y2D + x1Vn12D ↔ y2DVx1Vn12D
	(kon2, koff2), y2D + x2Vn21D ↔ y2DVx2Vn21D
	(kon2, koff2), y2D + x2Vn22D ↔ y2DVx2Vn22D
	(kon2, koff2), y1D + x1Vn11D ↔ y1DVx1Vn11D
	(kon2, koff2), y1D + x1Vn12D ↔ y1DVx1Vn12D
	(kon2, koff2), y1D + x2Vn21D ↔ y1DVx2Vn21D
	(kon2, koff2), y1D + x2Vn22D ↔ y1DVx2Vn22D
	(kon2, koff2), y2D + x1Vn11 ↔ y2DVx1Vn11
	(kon2, koff2), y2D + x2Vn21 ↔ y2DVx2Vn21
	(kon2, koff2), y2D + x1Vn12 ↔ y2DVx1Vn12
	(kon2, koff2), y2D + x2Vn22 ↔ y2DVx2Vn22
	(kon2, koff2), y1D + x1Vn11 ↔ y1DVx1Vn11
	(kon2, koff2), y1D + x2Vn21 ↔ y1DVx2Vn21
	(kon2, koff2), y1D + x1Vn12 ↔ y1DVx1Vn12
	(kon2, koff2), y1D + x2Vn22 ↔ y1DVx2Vn22
	(kon_p1, koff_p1), y1Vx1Vn11D + y2DVx1Vn11D ↔ y1Vx1Vn11D_y2DVx1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11D + y2DVx2Vn21D ↔ y1Vx1Vn11D_y2DVx2Vn21D
	(kon_p1, koff_p1), y1Vx1Vn11D + y1DVx1Vn11D ↔ y1Vx1Vn11D_y1DVx1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11D + y1DVx2Vn21D ↔ y1Vx1Vn11D_y1DVx2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21D + y2DVx1Vn11D ↔ y1Vx2Vn21D_y2DVx1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21D + y2DVx2Vn21D ↔ y1Vx2Vn21D_y2DVx2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21D + y1DVx1Vn11D ↔ y1Vx2Vn21D_y1DVx1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21D + y1DVx2Vn21D ↔ y1Vx2Vn21D_y1DVx2Vn21D
	(kon_p2, koff_p2), y2Vx1Vn12D + y2DVx1Vn12D ↔ y2Vx1Vn12D_y2DVx1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12D + y2DVx2Vn22D ↔ y2Vx1Vn12D_y2DVx2Vn22D
	(kon_p2, koff_p2), y2Vx1Vn12D + y1DVx1Vn12D ↔ y2Vx1Vn12D_y1DVx1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12D + y1DVx2Vn22D ↔ y2Vx1Vn12D_y1DVx2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22D + y2DVx1Vn12D ↔ y2Vx2Vn22D_y2DVx1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22D + y2DVx2Vn22D ↔ y2Vx2Vn22D_y2DVx2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22D + y1DVx1Vn12D ↔ y2Vx2Vn22D_y1DVx1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22D + y1DVx2Vn22D ↔ y2Vx2Vn22D_y1DVx2Vn22D
	kcat1, y1Vx1Vn11D_y2DVx1Vn11D → y1Vx1Vn11D + y2DVx1Vn11
	kcat1, y1Vx1Vn11D_y2DVx2Vn21D → y1Vx1Vn11D + y2DVx2Vn21
	kcat1, y1Vx1Vn11D_y1DVx1Vn11D → y1Vx1Vn11D + y1DVx1Vn11
	kcat1, y1Vx1Vn11D_y1DVx2Vn21D → y1Vx1Vn11D + y1DVx2Vn21
	kcat1, y1Vx2Vn21D_y2DVx1Vn11D → y1Vx2Vn21D + y2DVx1Vn11
	kcat1, y1Vx2Vn21D_y2DVx2Vn21D → y1Vx2Vn21D + y2DVx2Vn21
	kcat1, y1Vx2Vn21D_y1DVx1Vn11D → y1Vx2Vn21D + y1DVx1Vn11
	kcat1, y1Vx2Vn21D_y1DVx2Vn21D → y1Vx2Vn21D + y1DVx2Vn21
	kcat2, y2Vx1Vn12D_y2DVx1Vn12D → y2Vx1Vn12D + y2DVx1Vn12
	kcat2, y2Vx1Vn12D_y2DVx2Vn22D → y2Vx1Vn12D + y2DVx2Vn22
	kcat2, y2Vx1Vn12D_y1DVx1Vn12D → y2Vx1Vn12D + y1DVx1Vn12
	kcat2, y2Vx1Vn12D_y1DVx2Vn22D → y2Vx1Vn12D + y1DVx2Vn22
	kcat2, y2Vx2Vn22D_y2DVx1Vn12D → y2Vx2Vn22D + y2DVx1Vn12
	kcat2, y2Vx2Vn22D_y2DVx2Vn22D → y2Vx2Vn22D + y2DVx2Vn22
	kcat2, y2Vx2Vn22D_y1DVx1Vn12D → y2Vx2Vn22D + y1DVx1Vn12
	kcat2, y2Vx2Vn22D_y1DVx2Vn22D → y2Vx2Vn22D + y1DVx2Vn22
	(kon_p1, koff_p1), y1Vx1Vn11 + y2DVx1Vn11D ↔ y1Vx1Vn11_y2DVx1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11 + y2DVx2Vn21D ↔ y1Vx1Vn11_y2DVx2Vn21D
	(kon_p1, koff_p1), y1Vx1Vn11 + y1DVx1Vn11D ↔ y1Vx1Vn11_y1DVx1Vn11D
	(kon_p1, koff_p1), y1Vx1Vn11 + y1DVx2Vn21D ↔ y1Vx1Vn11_y1DVx2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21 + y2DVx1Vn11D ↔ y1Vx2Vn21_y2DVx1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21 + y2DVx2Vn21D ↔ y1Vx2Vn21_y2DVx2Vn21D
	(kon_p1, koff_p1), y1Vx2Vn21 + y1DVx1Vn11D ↔ y1Vx2Vn21_y1DVx1Vn11D
	(kon_p1, koff_p1), y1Vx2Vn21 + y1DVx2Vn21D ↔ y1Vx2Vn21_y1DVx2Vn21D
	(kon_p2, koff_p2), y2Vx1Vn12 + y2DVx1Vn12D ↔ y2Vx1Vn12_y2DVx1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12 + y2DVx2Vn22D ↔ y2Vx1Vn12_y2DVx2Vn22D
	(kon_p2, koff_p2), y2Vx1Vn12 + y1DVx1Vn12D ↔ y2Vx1Vn12_y1DVx1Vn12D
	(kon_p2, koff_p2), y2Vx1Vn12 + y1DVx2Vn22D ↔ y2Vx1Vn12_y1DVx2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22 + y2DVx1Vn12D ↔ y2Vx2Vn22_y2DVx1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22 + y2DVx2Vn22D ↔ y2Vx2Vn22_y2DVx2Vn22D
	(kon_p2, koff_p2), y2Vx2Vn22 + y1DVx1Vn12D ↔ y2Vx2Vn22_y1DVx1Vn12D
	(kon_p2, koff_p2), y2Vx2Vn22 + y1DVx2Vn22D ↔ y2Vx2Vn22_y1DVx2Vn22D
	kcat1, y1Vx1Vn11_y2DVx1Vn11D → y1Vx1Vn11 + y2DVx1Vn11
	kcat1, y1Vx1Vn11_y2DVx2Vn21D → y1Vx1Vn11 + y2DVx2Vn21
	kcat1, y1Vx1Vn11_y1DVx1Vn11D → y1Vx1Vn11 + y1DVx1Vn11
	kcat1, y1Vx1Vn11_y1DVx2Vn21D → y1Vx1Vn11 + y1DVx2Vn21
	kcat1, y1Vx2Vn21_y2DVx1Vn11D → y1Vx2Vn21 + y2DVx1Vn11
	kcat1, y1Vx2Vn21_y2DVx2Vn21D → y1Vx2Vn21 + y2DVx2Vn21
	kcat1, y1Vx2Vn21_y1DVx1Vn11D → y1Vx2Vn21 + y1DVx1Vn11
	kcat1, y1Vx2Vn21_y1DVx2Vn21D → y1Vx2Vn21 + y1DVx2Vn21
	kcat2, y2Vx1Vn12_y2DVx1Vn12D → y2Vx1Vn12 + y2DVx1Vn12
	kcat2, y2Vx1Vn12_y2DVx2Vn22D → y2Vx1Vn12 + y2DVx2Vn22
	kcat2, y2Vx1Vn12_y1DVx1Vn12D → y2Vx1Vn12 + y1DVx1Vn12
	kcat2, y2Vx1Vn12_y1DVx2Vn22D → y2Vx1Vn12 + y1DVx2Vn22
	kcat2, y2Vx2Vn22_y2DVx1Vn12D → y2Vx2Vn22 + y2DVx1Vn12
	kcat2, y2Vx2Vn22_y2DVx2Vn22D → y2Vx2Vn22 + y2DVx2Vn22
	kcat2, y2Vx2Vn22_y1DVx1Vn12D → y2Vx2Vn22 + y1DVx1Vn12
	kcat2, y2Vx2Vn22_y1DVx2Vn22D → y2Vx2Vn22 + y1DVx2Vn22
	syn_x1/(cell*avgd), ∅ ⇒ x1
	syn_x2/(cell*avgd), ∅ ⇒ x2
	syn_n11/(cell*avgd), ∅ ⇒ n11D
	syn_n12/(cell*avgd), ∅ ⇒ n12D
	syn_n21/(cell*avgd), ∅ ⇒ n21D
	syn_n22/(cell*avgd), ∅ ⇒ n22D
	syn_y1/(cell*avgd), ∅ ⇒ y1
	syn_y2/(cell*avgd), ∅ ⇒ y2
	deg_DHFR, x1Vn11D → ∅
	deg_DHFR, x1Vn12D → ∅
	deg_DHFR, x2Vn21D → ∅
	deg_DHFR, x2Vn22D → ∅
	deg_DHFR, y1Vx1Vn11D → ∅
	deg_DHFR, y1Vx1Vn12D → ∅
	deg_DHFR, y1Vx2Vn21D → ∅
	deg_DHFR, y1Vx2Vn22D → ∅
	deg_DHFR, y2Vx1Vn11D → ∅
	deg_DHFR, y2Vx1Vn12D → ∅
	deg_DHFR, y2Vx2Vn21D → ∅
	deg_DHFR, y2Vx2Vn22D → ∅
	deg_reg, n11 → ∅
	deg_reg, n21 → ∅
	deg_reg, n12 → ∅
	deg_reg, n22 → ∅
	deg_reg, x1Vn11 → ∅
	deg_reg, x2Vn21 → ∅
	deg_reg, x1Vn12 → ∅
	deg_reg, x2Vn22 → ∅
	deg_reg, y1Vx1Vn11 → ∅
	deg_reg, y1Vx2Vn21 → ∅
	deg_reg, y2Vx1Vn11 → ∅
	deg_reg, y2Vx2Vn21 → ∅
	deg_reg, y1Vx1Vn12 → ∅
	deg_reg, y1Vx2Vn22 → ∅
	deg_reg, y2Vx1Vn12 → ∅
	deg_reg, y2Vx2Vn22 → ∅
	deg_reg, y2D → ∅
	deg_reg, y1D → ∅
	deg_DHFR, y2DVx1Vn11D → ∅
	deg_DHFR, y2DVx1Vn12D → ∅
	deg_DHFR, y2DVx2Vn21D → ∅
	deg_DHFR, y2DVx2Vn22D → ∅
	deg_DHFR, y1DVx1Vn11D → ∅
	deg_DHFR, y1DVx1Vn12D → ∅
	deg_DHFR, y1DVx2Vn21D → ∅
	deg_DHFR, y1DVx2Vn22D → ∅
	deg_reg, y2DVx1Vn11 → ∅
	deg_reg, y2DVx2Vn21 → ∅
	deg_reg, y1DVx1Vn11 → ∅
	deg_reg, y1DVx2Vn21 → ∅
	deg_reg, y2DVx1Vn12 → ∅
	deg_reg, y2DVx2Vn22 → ∅
	deg_reg, y1DVx1Vn12 → ∅
	deg_reg, y1DVx2Vn22 → ∅
	deg_reg, x1 → ∅
	deg_reg, x2 → ∅
	deg_DHFR, n11D → ∅
	deg_DHFR, n12D → ∅
	deg_DHFR, n21D → ∅
	deg_DHFR, n22D → ∅
	deg_reg, y1 → ∅
	deg_reg, y2 → ∅
end kon1 koff1 kon2 koff2 kon_p1 koff_p1 kon_p2 koff_p2 kon_p3 koff_p3 deg_reg deg_DHFR  kcat1 kcat2 kcat3 cell avgd syn_x1 syn_x2 syn_n11 syn_n12 syn_n21 syn_n22 syn_y1 syn_y2;

p = [10^5,
100,
10^5,
10^-4,
10^5,
10^-4,
10^5,
10^-4,
10^5,
10^-4,
10^-5,
10^-3,
0.16,
0.16,
0.16,
4E-15,
6.02E23,
0.05,
0.06,
0.6,
0.8,
0.9,
0.6,
0.6,
0.6,
];

initial_condition = zeros(size(species(rn))[1])
tspan = (0., 1E6)
initial_condition_Gillespie = convert(Array{Int64}, initial_condition)
discrete_prob = DiscreteProblem(rn, initial_condition_Gillespie, tspan, p)
jump_prob = JumpProblem(rn, discrete_prob, RSSA())
@time sol = solve(jump_prob, SSAStepper(), saveat=10.)

[isaacsas]

I'm not sure why you are seeing no changes in the solution, however, I'm able to run to tspan=(0.0,1e12) just fine if I turn off saving the state after every jump:

jump_prob = JumpProblem(rn, discrete_prob, RSSA(), save_positions=(false,false))

(I also increased saveat to reduce memory use...)

[isaacsas]

Another problem I can see is with your zero order reactions. Somehow using the arrow that disables the mass action rate law is causing issues in the generated jumps. If you use a normal arrow like --> or → it looks like the right jumps get generated. (This still doesn't fix the nothing happening issue though.)

[isaacsas]

Are you sure you have your formula correct for the zero order reactions? The effective zero-order rates I'm seeing when I multiply out your parameters are on the order of 2e-10 to 2e-11, which means the average time till a reaction that actually makes anything would occur is ~10^10 to 10^11.

It might also be helpful with this many parameters to label them as follows:

params = Pair.(parameters(rn),p)

which gives

params = Pair.(parameters(rn),p)
25-element Vector{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Float64}}:
     kon1 => 100000.0
    koff1 => 100.0
     kon2 => 100000.0
    koff2 => 0.00010000000000000002
   kon_p1 => 100000.0
  koff_p1 => 0.00010000000000000002
   kon_p2 => 100000.0
  koff_p2 => 0.00010000000000000002
   kon_p3 => 100000.0
  koff_p3 => 0.00010000000000000002
  deg_reg => 1.0000000000000003e-5
 deg_DHFR => 0.0010000000000000002
    kcat1 => 0.16
    kcat2 => 0.16
    kcat3 => 0.16
     cell => 4.0e-15
     avgd => 6.02e23
   syn_x1 => 0.05
   syn_x2 => 0.06
  syn_n11 => 0.6
  syn_n12 => 0.8
  syn_n21 => 0.9
  syn_n22 => 0.6
   syn_y1 => 0.6
   syn_y2 => 0.6

This lets you check you've got them all setup correctly.

[uagaug]

Thanks for the trick on labeling parameters. These values should be correct as I was able to reproduce the deterministic ODE solution using

oprob = ODEProblem(rn, initial_condition, tspan, p)     
@time sol = solve(oprob, isoutofdomain=(u,p,t)->any(x->x<0,u), abstol=1e-12, reltol=1e-12)

which gives identical solutions compared to odeint in Python.

Your point about timescale is a great one, I am increasing tspan to 10^20 and will report back if it solves the issue.

[isaacsas]

Let's compare your parameters for x1. It is produced on a timescale of 10^10 units of time, but degrades on a timescale of 10^5. Just from that reaction we'd expect a steady-state amount of x1 of about 10^-10 / 10^-5 = 10^-5 for an ODE model, so in the stochastic model you'd expect to generally have none. Anyways, there could be a bug in Catalyst/DiffEqJump, but it seems like this might just be an issue of your chosen reaction rates. Are any of your ODE solutions actually larger than 1 when you look at them?

[uagaug]

Thank you. I followed this thought and played with the param values. It turns out the outputs from Gillespie now largely agrees with the outputs from deterministic ODE. It appears to be some scaling issues in my model. I will close this issue for now.

[isaacsas]

Glad we could help. Feel free to open another issue if you have other problems.