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

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.)

Open an issue in Catalyst.jl and we can follow up there. We’ll need to dive in and see what’s going on.

For posterity, the issue was Catalyst/DifferentialEqu and we decided that the problem was caused by using the arrows that disable rate laws for zero order reactions, saving state from all jumps, and with the choice of rate constants.