(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
(* Global constants                                                                               *)
(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
Na := 6.0221415e23 ;;
VolumeUnit := 1/Na ;;
fun gamma(V) = V ;; // actually Gamma(V)=V*VolumeUnit*Na, but VolumeUnit has been choosen in such a way that Gamma(V)=V waiting for the right volume (mean volume of a germinal bacterium)

K_CreD   := 0.2   ;;  // Cre degradation
K_DAPiD  := 0.2   ;;  // DAP degradation
K_DAPApI := 0.5   ;;  // DAPAp inhibition by DAP
K_DAPApA := 0.1   ;;  // DAPAp desinhibition releasing a DAP molecule
K_CreP1  := 1.0   ;;  // Cre production (DAPAp not inhibited)
K_CreP2  := 0.001 ;;  // Cre production (DAPAp inhibited)
K_Diff   := 0.05  ;;  // Differentiation
K_DAPiP  := 1.0   ;;  // DAP production
K_DAPEx  := 0.25   ;;  // DAP export
K_DAPIm  := 0.5   ;;  // DAP import
K_DAPeD  := 0.1   ;;  // DAP degradation

K_Div    := 0.02   ;; // Mitosis rate
K_Death  := 0.01   ;; // Death rate

fun fstOrder(k,v) = k ;;
fun sndOrder(k,v) = k/gamma(v) ;;

appliedBactRule := `Unkown  ;;
DAPEx := 0;
trans inBact[DAPE,v=1] = {

  (* Chemicals degradation *)
  `Cre ={ C = fstOrder(K_CreD,1.0)  }=> (appliedBactRule := "CreD"  ; <undef>) ;
  `DAP ={ C = fstOrder(K_DAPiD,1.0) }=> (appliedBactRule := "DAPiD" ; <undef>) ;

  (* Promotor activity *)
  `DAPAp, `DAP ={ C = sndOrder(K_DAPApI,1.0) }=> (appliedBactRule := "DAPApI" ; `DAPAp_i) ;
  `DAPAp_i     ={ C = fstOrder(K_DAPApA,1.0) }=> (appliedBactRule := "DAPApA" ; `DAPAp, `DAP) ;
  `DAPAp       ={ C = fstOrder(K_CreP1,1.0)  }=> (appliedBactRule := "CreP1"  ; `DAPAp, `Cre) ;
  `DAPAp_i     ={ C = fstOrder(K_CreP2,1.0)  }=> (appliedBactRule := "CreP2"  ; `DAPAp_i, `Cre) ;

  (* Differentiation (exclusively for germinal cells) *)
  `LOXP_Box, `Cre     ={ C = sndOrder(2.0*K_Diff,1.0) }=> (appliedBactRule := "Diff1" ; `LOXP_Box_Cre) ;
  `LOXP_Box_Cre, `Cre ={ C = sndOrder(K_Diff,1.0)     }=> (appliedBactRule := "Diff2" ; `DAP_Box) ;

  (* DAP production (exclusively for somatic cells) *)
  `DAP_Box ={ C = fstOrder(K_DAPiP,1.0) }=> (appliedBactRule := "DAPiP" ; `DAP_Box, `DAP) ;

  (* DAP externalization & internalization *)
  `DAP   ={ C = fstOrder(K_DAPEx,1.0)                }=> (appliedBactRule := "DAPEx"; DAPEx := DAPEx+v; <undef>) ;
  _ as x ={ A = (\c.(DAPEx * fstOrder(K_DAPIm,1.0))) }=> (appliedBactRule := "DAPIm"; DAPEx := DAPEx-v; x, `DAP) ;

} ;;

output := "/tmp/smb.dat" ;;
fun SMBoutput[ofile=output](m,c) = (
  if (iteration%m == 0) then (
    stdout << iteration << " " << tau << "\n";
    ofile           <<
    tau             << " " << 
    appliedBactRule << " " << 
    DAPEx           << " " << 
    count(`DAP,c)   << " " << 
    count(`Cre,c)   << "\n"
  ) else 0 fi;
  c
) ;;

DAPEx:=1000;
inBact[DAPE=10,
       v=0,
       strategy=`gillespie,
       postlude=SMBoutput(10),
       interlude=SMBoutput(10),
       prelude=SMBoutput(10),
       fixpoint=(\c1.\c2.(tau>=100))
      ](bagify(((#10 `DAP), `DAPAp, `LOXP_Box, seq:())));;

system("gnuplot -persist script") ;;