Circuit with program regulation

circuit_with_program_reg.ka

// Model to run by KappaSim class
// 
// Agents
//
// A and B can be phosphorylated/unphosphorylated
%agent: A(head, tail)
%agent: B(head, tail)

// head of transition counter polymer 
%agent: T(tail)

// memory molecules
%agent: AToA(site)
%agent: AToB(site)
%agent: BToA(site)
%agent: BToB(site)

// memory molecule-associated integrators
%agent: AToBIntegrator(site, nutr_site)
%agent: BToAIntegrator(site, nutr_site)

// gene expression programs for each state
%agent: AProgram()
%agent: BProgram()

//
// Rules
//

// free A or B can bind (irreversibly) to transition polymer starting
// monomer which I will call T
//  (this is the 'base case', so to speak, of polymerization)
//new rule
T(tail[.]),A(head[.],tail[.]) -> T(tail[1]),A(head[1],tail[.]) @ 'nutr_bind_polystart_k'
T(tail[.]),B(head[.],tail[.]) -> T(tail[1]),B(head[1],tail[.]) @ 'nutr_bind_polystart_k'

// the 'inductive step' so to speak of polymerization: A or B with a free head
// can bind to an A or B (irreversibly) which has free tail but a *bound* head 
// (which can only happen if already bound to T)
// when this happens, the appropriate transition memory molecules
// are produced ex nihilo.
// four cases here: 
// (1) A -> A
A(head[_],tail[.]),A(head[.],tail[.]),. -> A(head[_],tail[1]),A(head[1],tail[.]),AToA() @ 'nutr_add_to_poly_k'
// (2) A -> B 
A(head[_],tail[.]),B(head[.],tail[.]),. -> A(head[_],tail[1]),B(head[1],tail[.]),AToB() @ 'nutr_add_to_poly_k'
// (3) B -> A
B(head[_],tail[.]),A(head[.],tail[.]),. -> B(head[_],tail[1]),A(head[1],tail[.]),BToA() @ 'nutr_add_to_poly_k'
// (4) B -> B
B(head[_],tail[.]),B(head[.],tail[.]),. -> B(head[_],tail[1]),B(head[1],tail[.]),BToB() @ 'nutr_add_to_poly_k'

//// AToB and AToA (reversibly) bind, in mutually exclusive fashion,
//// to AToBIntegrator
//// ('communist binding')
// note that AToB affinity to Integrator needs to be X times greater than AToB's 
// to compensate for fact that we only add one B to the polymer for an AToB transition
// but we add multiple As when we are in a stretch of As... so the X-factor depends
// on the duration (and size) of the nutrient pulse
AToB(site[.]),AToBIntegrator(site[.]) <-> AToB(site[1]),AToBIntegrator(site[1]) @ 'AToB_bind_int_for_k', 'AToB_bind_int_rev_k'
AToA(site[.]),AToBIntegrator(site[.]) <-> AToA(site[1]),AToBIntegrator(site[1]) @ 'AToA_bind_int_for_k', 'AToA_bind_int_rev_k'

//// BToA and BToB (reversibly) bind, in mutually exclusive fashion,
//// to BToAIntegrator
//// ('communist binding')
// note that BToA affinity to Integrator needs to be X times greater than BToA's 
// to compensate for fact that we only add one A to the polymer for an BToA transition
// but we add multiple Bs when we are in a stretch of Bs... so the X-factor depends
// on the duration (and size) of the nutrient pulse
BToA(site[.]),BToAIntegrator(site[.]) <-> BToA(site[1]),BToAIntegrator(site[1]) @ 'BToA_bind_int_for_k', 'BToA_bind_int_rev_k'
BToB(site[.]),BToAIntegrator(site[.]) <-> BToB(site[1]),BToAIntegrator(site[1]) @ 'BToB_bind_int_for_k', 'BToB_bind_int_rev_k'

// integrator bound to memory molecule can bind the tail of a free nutrient
AToBIntegrator(site[_],nutr_site[.]),A(head[.],tail[.]) -> AToBIntegrator(site[_],nutr_site[1]),A(head[.],tail[1]) @ 'int_bind_nutr_for_k'
// integrator and nutrient can dissociate (regardless of whether integrator is bound to memory molecule)
AToBIntegrator(nutr_site[1]),A(head[.],tail[1]) -> AToBIntegrator(nutr_site[.]),A(head[.],tail[.]) @ 'int_bind_nutr_rev_k'

///
/// regulation of AToB integrator
///
/*
* rules for active AToB integrator complex 
*/
// active complex (AToB,AToBIntegrator,A) produces BProgram 
AToB(site[1]),AToBIntegrator(site[1],nutr_site[2]),A(head[.],tail[2]),. -> AToB(site[1]),AToBIntegrator(site[1],nutr_site[2]),A(head[.],tail[2]),BProgram() @ 'ActAToBInt_BProgram_on_k'

/*
* rules for repressed AToA integrator complex 
*/
// repressed complex (AToA,AToBIntegrator,A) degrades BProgram 
AToA(site[1]),AToBIntegrator(site[1],nutr_site[2]),A(head[.],tail[2]),BProgram() -> AToA(site[1]),AToBIntegrator(site[1],nutr_site[2]),A(head[.],tail[2]),. @ 'RepAToBInt_BProgram_off_k'

// integrator bound to memory molecule can bind the tail of a free nutrient
BToAIntegrator(site[_],nutr_site[.]),B(head[.],tail[.]) -> BToAIntegrator(site[_],nutr_site[1]),B(head[.],tail[1]) @ 'int_bind_nutr_for_k'
// integrator and nutrient can dissociate (regardless of whether integrator is bound to memory molecule)
BToAIntegrator(nutr_site[1]),B(head[.],tail[1]) -> BToAIntegrator(nutr_site[.]),B(head[.],tail[.]) @ 'int_bind_nutr_rev_k'

///
/// regulation of BToA integrator
///
/*
* rules for active BToA integrator complex 
*/
// active complex (BToA,BToAIntegrator,B) degrades BProgram 
BToA(site[1]),BToAIntegrator(site[1],nutr_site[2]),B(head[.],tail[2]),BProgram() -> BToA(site[1]),BToAIntegrator(site[1],nutr_site[2]),B(head[.],tail[2]),. @ 'ActBToAInt_BProgram_off_k'

/*
* rules for repressed BToB integrator complex 
*/
// repressed complex (BToB,BToAIntegrator,B) produces BProgram 
BToB(site[1]),BToAIntegrator(site[1],nutr_site[2]),B(head[.],tail[2]),. -> BToB(site[1]),BToAIntegrator(site[1],nutr_site[2]),B(head[.],tail[2]),BProgram() @ 'RepBToAInt_BProgram_on_k'

/* Degradation reactions */

// degrade nutrients (only unbound nutrients can be degraded)
'A_deg' A(head[.],tail[.]) -> . @ 'nutr_deg_k'

// degrade gene expression programs
'AProgram_deg' AProgram() -> . @ 'program_deg_k'

// degrade nutrients (only unbound nutrients can be degraded)
'B_deg' B(head[.],tail[.]) -> . @ 'nutr_deg_k'

// degrade gene expression programs
'BProgram_deg' BProgram() -> . @ 'program_deg_k'

// degrade memory molecules
'AToB_deg' AToB() -> . @ 'mem_deg_k'
'AToA_deg' AToA() -> . @ 'mem_deg_k'

// degrade memory molecules
'BToA_deg' BToA() -> . @ 'mem_deg_k'
'BToB_deg' BToB() -> . @ 'mem_deg_k'

%var: 'nutr_bind_polystart_k' 0.0001
%var: 'nutr_add_to_poly_k' 0.0001
%var: 'AToB_bind_int_for_k' 0.000337746096401907
%var: 'AToB_bind_int_rev_k' 0.005
%var: 'AToA_bind_int_for_k' 0.0000183457482115204
%var: 'AToA_bind_int_rev_k' 0.005
%var: 'BToA_bind_int_for_k' 0.000337746096401907
%var: 'BToA_bind_int_rev_k' 0.005
%var: 'BToB_bind_int_for_k' 0.0000183457482115204
%var: 'BToB_bind_int_rev_k' 0.005
%var: 'nutr_deg_k' 1e-05
%var: 'mem_deg_k' 1e-05
%var: 'program_deg_k' 0.0
%var: 'int_bind_nutr_for_k' 0.005
%var: 'int_bind_nutr_rev_k' 0.005
%var: 'ActAToBInt_BProgram_on_k' 0.001
%var: 'RepAToBInt_BProgram_off_k' 0.001
%var: 'ActBToAInt_BProgram_off_k' 0.001
%var: 'RepBToAInt_BProgram_on_k' 0.001
%var: 'RepAToBInt_AProgram_on_k' 'ActAToBInt_BProgram_on_k'
%var: 'ActAToBInt_AProgram_off_k' 'RepAToBInt_BProgram_off_k'
%var: 'ActBToAInt_AProgram_on_k' 'RepBToAInt_BProgram_on_k'
%var: 'RepBToAInt_AProgram_off_k' 'ActBToAInt_BProgram_off_k'
%var: 'amount_int' 400
%var: 'amount_T' 100
%var: 'amount_B' 0
%var: 'amount_A' 0
%var: 'amount_AToA' 0
%var: 'amount_AToB' 0
%var: 'amount_BToB' 0
%var: 'amount_BToA' 0
%init: amount_int AToBIntegrator()
%init: amount_AToA AToA()
%init: amount_AToB AToB()
%init: amount_int BToAIntegrator()
%init: amount_BToB BToB()
%init: amount_BToA BToA()
%init: amount_B B()
%init: amount_A A()
%init: amount_T T()

%obs: 'A' |A(head[.], tail[.])|
%obs: 'B' |B(head[.], tail[.])|
%obs: 'AToA' |AToA()|
%obs: 'AToB' |AToB()|
%obs: 'BToA' |BToA()|
%obs: 'BToB' |BToB()|
%obs: 'AProgram' |AProgram()|
%obs: 'BProgram' |BProgram()|
%obs: 'ActiveAToBIntegrator' |AToBIntegrator(site[1]),AToB(site[1])|
%obs: 'RepressedAToBIntegrator' |AToBIntegrator(site[1]),AToA(site[1])|
%obs: 'ActiveBToAIntegrator' |BToAIntegrator(site[1]),BToA(site[1])|
%obs: 'RepressedBToAIntegrator' |BToAIntegrator(site[1]),BToB(site[1])|
%obs: 'free_AToBIntegrator' |AToBIntegrator(site[.])|
%obs: 'free_BToAIntegrator' |BToAIntegrator(site[.])|
%obs: 'AToBIntegrator' |AToBIntegrator()|
%obs: 'BToAIntegrator' |BToAIntegrator()|
%obs: 'bound_AToBIntegrator' |AToBIntegrator(site[_])|
%obs: 'bound_BToAIntegrator' |BToAIntegrator(site[_])|
%obs: 'IntAToB' |AToBIntegrator(site[1]),AToB(site[1])|
%obs: 'IntAToA' |AToBIntegrator(site[1]),AToA(site[1])|
%obs: 'IntBToA' |BToAIntegrator(site[1]),BToA(site[1])|
%obs: 'IntBToB' |BToAIntegrator(site[1]),BToB(site[1])|
%obs: 'A_IntAToB' |AToBIntegrator(site[1],nutr_site[2]),AToB(site[1]),A(tail[2])|
%obs: 'A_IntAToA' |AToBIntegrator(site[1],nutr_site[2]),AToA(site[1]),A(tail[2])|
%obs: 'B_IntBToA' |BToAIntegrator(site[1],nutr_site[2]),BToA(site[1]),B(tail[2])|
%obs: 'B_IntBToB' |BToAIntegrator(site[1],nutr_site[2]),BToB(site[1]),B(tail[2])|

        // pulse size for A
        %var: 'A_pulse_size' 5000
        
        // pulse size for B
        %var: 'B_pulse_size' 5000
        %mod: [true] do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: [true] do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 49.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 50.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 99.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 100.0 do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 149.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 150.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 199.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 200.0 do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 249.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 250.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 299.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 300.0 do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 349.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 350.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 399.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 400.0 do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 449.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 450.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 499.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 500.0 do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 549.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 550.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 599.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 600.0 do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 649.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 650.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 699.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 700.0 do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 749.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 750.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 799.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 800.0 do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 849.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 850.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 899.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 900.0 do $ADD 'A_pulse_size' A(head[.],tail[.]);
%mod: alarm 949.999 do $DEL |A(head[.],tail[.])| A(head[.],tail[.]);$DEL |B(head[.],tail[.])| B(head[.],tail[.]);%mod: alarm 950.0 do $ADD 'B_pulse_size' B(head[.],tail[.]);
%mod: alarm 999.0 do $SNAPSHOT "snap".[E].".json";

%mod: ([T] > 1000.0) do $STOP;

// perturbations