directive duration 7000.0 points 1000
directive plot sum(<_ X1^ X2 X3^>); sum(<_ Y1^ Y2 Y3^>); sum(<_ Z1^ Z2 Z3^>)

(* First define macros for the CRN->DNA translation scheme *)
(* creates a DNA signal species for CRN species A *)
def signal(N,a1,a2,a3) = new h (N * <h a1^ a2 a3^>)

(* creates DNA fuels for a unimolecular reaction A -> B *)
def reaction1to1(N,a1,a2,a3,b1,b2,b3) =
  new r1
  ( N * {a1^*}[a2 a3^]<r1 b1^>
  | N * {a3^*}[r1 b1^]<b2 b3^> )

(* creates DNA fuels for a unimolecular reaction A -> B+C *)
def reaction1to2(N,a1,a2,a3,b1,b2,b3,c1,c2,c3) =
  new r1
  ( N * {a1^*}[a2 a3^]<r1 b1^ r2 c1^>
  | N * {a3^*}[r1 b1^]<b2 b3^>:[r2 c1^]<c2 c3^> )

(* creates DNA fuels for a bimolecular reaction A+B -> C *)
def reaction2to1(N,a1,a2,a3,b1,b2,b3,c1,c2,c3) =
  new r1
  ( N * {a1^*}[a2 a3^ b1^]:[b2 b3^]<r1 c1^>
  | N * <a2 a3^ b1^>
  | N * {b3^*}[r1 c1^]<c2 c3^> )

(* creates DNA fuels for a bimolecular reaction A+B -> C+D *)
def reaction2to2(N,a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3) =
  new r1 new r2
  ( N * {a1^*}[a2 a3^ b1^]:[b2 b3^]<r1 c1^ r2 d1^>
  | N * <a2 a3^ b1^>
  | N * {b3^*}[r1 c1^]<c2 c3^>:[r2 d1^]<d2 d3^> )

(* Now populate the test tube with DNA fuels and signal strands *)
( reaction2to2(1000, X1,X2,X3, Y1,Y2,Y3, Y1,Y2,Y3, Y1,Y2,Y3)
| reaction2to2(1000, Y1,Y2,Y3, Z1,Z2,Z3, Z1,Z2,Z3, Z1,Z2,Z3)
| reaction2to2(1000,  Z1,Z2,Z3, X1,X2,X3, X1,X2,X3, X1,X2,X3)
| signal(200, X1,X2,X3)
| signal(40, Y1,Y2,Y3)
| signal(20, Z1,Z2,Z3)
)