Paris/Sources

From 2007.igem.org

(Difference between revisions)
(→Global Variables)
(→Structure)
Line 435: Line 435:
==Structure==
==Structure==
 +
<code><pre>
 +
(* ********************************************************************************************** *)
 +
(* ********************************************************************************************** *)
 +
(* State of the system                                                                            *)
 +
(* ********************************************************************************************** *)
 +
(* ********************************************************************************************** *)
 +
(* The system is composed of bacteria of different kinds: germinal and somatic                    *)
 +
(*                                                                                                *)
 +
(* It is modeled by nested multi-sets. The super multi-set corresponds to the environment where  *)
 +
(* chemicals and bacteria are diffusing.                                                          *)
 +
(* A set of functions is added to memorize the quantity of chemicals in the system in order not  *)
 +
(* to evaluate it at each time.                                                                  *)
 +
(* ********************************************************************************************** *)
 +
(* ********************************************************************************************** *)
 +
 +
cptCre  := 0 ;;
 +
cptDAPe  := 0 ;;
 +
cptDAPi  := 0.0 ;;
 +
cptDAP  := 0 ;;
 +
cptBactS := 0 ;;
 +
cptBactG := 0 ;;
 +
cptBact  := 0 ;;
 +
fun incrCre() = (cptCre := cptCre + 1) ;;
 +
fun decrCre() = (cptCre := cptCre - 1) ;;
 +
fun incrDAP() = (cptDAP := cptDAP + 1) ;;
 +
fun decrDAP() = (cptDAP := cptDAP - 1) ;;
 +
fun incrDAPe() = (cptDAPe := cptDAPe + 1) ;;
 +
fun decrDAPe() = (cptDAPe := cptDAPe - 1) ;;
 +
fun incrDAPi() = (cptDAPi := cptDAPi + 1) ;;
 +
fun decrDAPi() = (cptDAPi := cptDAPi - 1) ;;
 +
fun DAPe_DAPi() = (incrDAPi() ; decrDAPe()) ;;
 +
fun DAPi_DAPe() = (incrDAPe() ; decrDAPi()) ;;
 +
fun incrNCre(N) = (cptCre := cptCre + N) ;;
 +
fun decrNCre(N) = (cptCre := cptCre - N) ;;
 +
fun incrNDAP(N) = (cptDAP := cptDAP + N) ;;
 +
fun decrNDAP(N) = (cptDAP := cptDAP - N) ;;
 +
fun incrNDAPe(N) = (cptDAPe := cptDAPe + N) ;;
 +
fun decrNDAPe(N) = (cptDAPe := cptDAPe - N) ;;
 +
fun incrNDAPi(N) = (cptDAPi := cptDAPi + N) ;;
 +
fun decrNDAPi(N) = (cptDAPi := cptDAPi - N) ;;
 +
 +
fun incrBact() = (cptBact := cptBact + 1) ;;
 +
fun decrBact() = (cptBact := cptBact - 1) ;;
 +
fun incrBactS() = (cptBactS := cptBactS + 1) ;;
 +
fun decrBactS() = (cptBactS := cptBactS - 1) ;;
 +
fun incrBactG() = (cptBactG := cptBactG + 1) ;;
 +
fun decrBactG() = (cptBactG := cptBactG - 1) ;;
 +
fun BactS_BactG() = (raise  (`Erreur "Impossible S to G differentiation.")) ;;
 +
fun BactG_BactS() = (incrBactS() ; decrBactG()) ;;
 +
fun incrNBact(N) = (cptBact := cptBact + N) ;;
 +
fun decrNBact(N) = (cptBact := cptBact - N) ;;
 +
fun incrNBactS(N) = (cptBactS := cptBactS + N) ;;
 +
fun decrNBactS(N) = (cptBactS := cptBactS - N) ;;
 +
fun incrNBactG(N) = (cptBactG := cptBactG + N) ;;
 +
fun decrNBactG(N) = (cptBactG := cptBactG - N) ;;
 +
 +
collection Bact      = bag
 +
and constraint BactG = [~`DAP_Box]Bact
 +
and constraint BactS = [~`LOXP_Box && ~`LOXP_Box_Cre]Bact
 +
and collection Env  = bag  ;;
 +
 +
id := 0 ;;
 +
fun new_bactId() = (id := id+1 ; id) ;;
 +
fun new_bact(t)  = (incrBact() ; Bactify({birth=t, date=0.0, id=new_bactId(), next_state=`Unknown}::seq:()) )  ;;
 +
fun new_bactS(t) = (incrBactS() ; `DAPAp :: `DAP_Box  :: new_bact(t)) ;;
 +
fun new_bactG(t) = (incrBactG() ; (*incrDAPi() ; incrDAP() ;*) `DAPAp :: `LOXP_Box :: new_bact(t)) ;;
 +
fun new_env(t, nbG, nbS) = (
 +
  incrNDAP(nb_dape) ;
 +
  incrNDAPe(nb_dape) ;
 +
  let env = fold((\i.\acc.((new_bactG(t))::acc)),Envify(fold((\n.\acc.(`DAP::acc)),seq:(),nb_dape)), nbG) in
 +
    fold((\i.\acc.((new_bactS(t))::acc)),env, nbS)
 +
) ;;
 +
fun map_bact(f,env) = (
 +
  map((\e.(if Bact(e) then f(e) else e fi)),env) //update_bact
 +
) ;;
 +
 +
fun divide_bactG(t,b,id) = (
 +
  let id' = new_bactId() in
 +
  let nbDAP = count(`DAP,b)
 +
  and nbCre = count(`Cre,b) in
 +
  let bact1_core = Bactify({birth=t, date=0.0, id=id,  next_state=`Unknown}::seq:())
 +
  and bact2_core = Bactify({birth=t, date=0.0, id=id', next_state=`Unknown}::seq:()) in
 +
  let bact1_core = fold((\n.\acc.(`DAP::acc)),bact1_core,nbDAP/2)
 +
  and bact2_core = fold((\n.\acc.(`DAP::acc)),bact2_core,nbDAP - nbDAP/2) in
 +
  let bact2_core = fold((\n.\acc.(`Cre::acc)),bact2_core,nbCre/2)
 +
  and bact1_core = fold((\n.\acc.(`Cre::acc)),bact1_core,nbCre - nbCre/2) in
 +
  let bact1_core = if (member(`DAPAp_i,b)) then (`DAPAp_i) else (`DAPAp) fi :: bact1_core
 +
  and bact2_core = `DAPAp :: bact2_core in
 +
  let bact2_core = if (member(`LOXP_Box_Cre,b)) then (`LOXP_Box_Cre) else (`LOXP_Box) fi :: bact2_core
 +
  and bact1_core = `LOXP_Box :: bact1_core in
 +
    bact1_core :: bact2_core :: seq:()
 +
) ;;
 +
 +
fun set_DAPe(N,env) = (
 +
  let dDAP = cptDAPe - N in
 +
    if (dDAP > 0)
 +
    then (
 +
      (* cptDAPe > N : trop de DAP *)
 +
      cptDAPe := N ;
 +
      cptDAP  := cptDAP - dDAP ;
 +
      diff(env, fold((\n.\acc.(`DAP::acc)),Env:(),dDAP))
 +
    )
 +
    else (
 +
      if (dDAP < 0)
 +
      then (
 +
(* cptDAPe < N : pas assez de DAP *)
 +
cptDAPe := N ;
 +
cptDAP  := cptDAP - dDAP ;
 +
join(env, fold((\n.\acc.(`DAP::acc)),Env:(),-1*dDAP))
 +
      )
 +
      else (env) fi
 +
    ) fi
 +
) ;;
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
fun static_help() = (
 +
  stdout << "\n** Help of the \"static\" module **\n\n" ;
 +
  stdout << "Globals\n" ;
 +
  stdout << "\t- cptCre (cptCre): `Cre counter\n\n" ;
 +
  stdout << "\t- cptDAP (cptDAP): `DAP counter\n\n" ;
 +
  stdout << "\t- cptDAPe (cptDAPe): `DAP counter (in the environment)\n\n" ;
 +
  stdout << "\t- cptDAPi (cptDAPi): `DAP counter (in bacteria)\n\n" ;
 +
  stdout << "\t- cptBact (cptBact): bacteria counter\n\n" ;
 +
  stdout << "\t- cptBactG (cptBactG): germinal bacteria counter\n\n" ;
 +
  stdout << "\t- cptBactS (cptBactS): somatic bacteria counter\n\n" ;
 +
 +
  stdout << "Functions\n" ;
 +
  stdout << "\t- incrCre() (unit -> int):\n\t\tincrement cptCre\n\n" ;
 +
  stdout << "\t- decrCre() (unit -> int):\n\t\tdecrement cptCre\n\n" ;
 +
  stdout << "\t- incrDAP() (unit -> int):\n\t\tincrement cptDAP\n\n" ;
 +
  stdout << "\t- decrDAP() (unit -> int):\n\t\tdecrement cptDAP\n\n" ;
 +
  stdout << "\t- incrDAPe() (unit -> int):\n\t\tincrement cptDAPe\n\n" ;
 +
  stdout << "\t- decrDAPe() (unit -> int):\n\t\tdecrement cptDAPe\n\n" ;
 +
  stdout << "\t- incrDAPi() (unit -> int):\n\t\tincrement cptDAPi\n\n" ;
 +
  stdout << "\t- decrDAPi() (unit -> int):\n\t\tdecrement cptDAPi\n\n" ;
 +
  stdout << "\t- DAPe_DAPi() (unit -> int):\n\t\tdecrement cptDAPe, increment cptDAPi\n\n" ;
 +
  stdout << "\t- DAPi_DAPe() (unit -> int):\n\t\tdecrement cptDAPi, increment cptDAPe\n\n" ;
 +
  stdout << "\t- incrBact() (unit -> int):\n\t\tincrement cptBact\n\n" ;
 +
  stdout << "\t- decrBact() (unit -> int):\n\t\tdecrement cptBact\n\n" ;
 +
  stdout << "\t- incrBactG() (unit -> int):\n\t\tincrement cptBactG\n\n" ;
 +
  stdout << "\t- decrBactG() (unit -> int):\n\t\tdecrement cptBactG\n\n" ;
 +
  stdout << "\t- incrBactS() (unit -> int):\n\t\tincrement cptBactS\n\n" ;
 +
  stdout << "\t- decrBactS() (unit -> int):\n\t\tdecrement cptBactS\n\n" ;
 +
  stdout << "\t- BactG_BactS() (unit -> int):\n\t\tdecrement cptBactG, increment cptBactS\n\n" ;
 +
 +
  stdout << "\t- new_Bact(t) (float -> Bact):\n\t\tcreate a new bacterium at time t\n\n" ;
 +
  stdout << "\t- new_BactG(t) (float -> BactG):\n\t\tcreate a new germinal bacterium at time t\n\n" ;
 +
  stdout << "\t- new_BactS(t) (float -> BactS):\n\t\tcreate a new somatic bacterium at time t\n\n" ;
 +
  stdout << "\t- new_env(t,ng,ns) (float * int * int -> Env):\n\t\tcreate a new environment containing ng germinal bacteria and ns somatic bacteria at time t\n\n" ;
 +
  stdout << "\t- map_bact(f,env) ((Bact -> Bact) * Env -> Env):\n\t\tapply the function f on each bacterium of env\n\n" ;
 +
  stdout << "\t- set_DAPe(N,env) (int * Env -> Env):\n\t\tchange env, cptDAP and cptDAPe to N `DAP\n\n"
 +
) ;;
 +
</pre></code>
==Scheduling==
==Scheduling==

Revision as of 18:21, 24 October 2007



Contents

Cell auto

!include "gbviewOutput.mgs" ;;

N := 30;;

gbf grille = <nord, est ; 30 * nord, 30 * est>
and record Bact = { DAP:float, DAPe:float }
and record BactG = Bact + { bactG }
and record BactS = Bact + { bactS} ;;

fun iota_v(v,n) = map((\x.v), iota(n,seq:())) ;;
fun new_grid(v) = iota_v(iota_v(v,N), N) following |nord>, |est> ;;

fun majBactG[DAPexport = 0.2, DAPimport = 0.2, DAPconso= 1, SupDiff = 2.4, InfDiff = 2.02 ](x) = 
 	if ( ( x.DAP > InfDiff ) & ( x.DAP < SupDiff ) & (random(1) < 0.8) ) then
		{DAP = x.DAP , DAPe = x.DAPe, bactS} else
		let dap = if x.DAP - DAPconso  > 0 then  x.DAP - DAPconso  else  0 fi
		 and dape= x.DAPe
		 in if 	dap > dape then 
                    x + {DAP = dap - DAPexport * dap, DAPe = dape + DAPexport * dap }
		   else x + {DAP = dap + DAPimport * dape, DAPe = dape - DAPimport * dape }
                  fi
        fi 
;;

fun majBactS[DAPexport = 0.2, DAPimport = 0.2, DAPprod=10](x) =  
		let dap =  x.DAP + DAPprod
		 and dape= x.DAPe
		 in if 	dap > dape then 
                    x + {DAP = dap - DAPexport * dap, DAPe = dape + DAPexport * dap }
		   else x + {DAP = dap + DAPimport * dape, DAPe = dape - DAPimport * dape }
                  fi

;;


trans evol[Delta_t=0.1,DAPeDiff=1, DAPeDegrad = 0.2] = {



  (* production of DAP *)
  x:BactS => ( 
	let d = neighborsfold(
      (\y.\acc.( DAPeDiff*Delta_t*(y.DAPe-x.DAPe) + acc)),
      x.DAPe,
      x)
    in majBactS( x + {DAPe = if ( d - DAPeDegrad * x.DAPe) > 0 then d - DAPeDegrad * x.DAPe else 0 fi } )
      );

  (* Diffirentiation and Consumption *)
 x:BactG => ( 
	let d = neighborsfold(
      (\y.\acc.( DAPeDiff*Delta_t*(y.DAPe-x.DAPe) + acc)),
      x.DAPe,
      x)
    in majBactG(x + {DAPe = if ( d - DAPeDegrad * x.DAPe) > 0 then  d - DAPeDegrad * x.DAPe else 0 fi } )
    

  )

} ;;





g := new_grid({DAP=6.0, DAPe=3.0, bactG}) ; 0 ;;
g := set_gbfpos(g, (random(30)*|nord> + random(30)*|est>), {DAP = 6.0, DAPe=3.0, bactS}) ;;
g := set_gbfpos(g, (random(30)*|nord> + random(30)*|est>), {DAP = 6.0, DAPe=3.0, bactS}) ;;
g := set_gbfpos(g, (random(30)*|nord> + random(30)*|est>), {DAP = 6.0, DAPe=3.0, bactS}) ;;
g := set_gbfpos(g, (random(30)*|nord> + random(30)*|est>), {DAP = 6.0, DAPe=3.0, bactS}) ;; 
evol[iter=1000,interlude=GBVexport((\v.("DAPe = "+(12 * v.DAPe )+if BactS(v) then ", bactS= 255" else  ", bactG = 255"fi )))](g) ;;



Cell auto 2

Constant rate of differentiation

record MecaBact = {x, y, vx, vy, fx, fy, radius}
and record CellBact  = {dap:float, soma:bool}
and record Bact = MecaBact + CellBact;;


delaunay(2) D2 = (\e.(e.x, e.y)) ;;

fun noise(x) = x + 0.0005 - random(0.001) ;;
// --- Mechanistic-------------------------------------------------------

DT    := 0.05;;
K     := 1.0;;
MU    := 1.8;;
R0_Gm := 0.50;; //germinal cells minmal rayon 
R0_G  :=0.75;;
R0_S  := 1.00;;

// interaction : the effect of the cells between each others
fun interaction(ref, src) = (
 let X      = ref.x - src.x
 and Y      = ref.y - src.y in
 let dist   = sqrt(X*X+Y*Y) in
 let spring = 0.0-K*(dist-(ref.radius+src.radius))/dist in
   {fx=X*spring - ref.vx*MU, fy = Y*spring - ref.vy*MU}
) ;;

fun add_vect(u, v) = { fx = u.fx + v.fx, fy = u.fy + v.fy } ;;
fun sum(x, u, acc) = add_vect(acc, interaction(x,u)) ;;


trans Meca = {
 e => (
   let f = neighborsfold(sum(e), {fx=0,fy=0}, e) in
     e+{ x  = noise(e.x  + DT*e.vx),
         y  = noise(e.y  + DT*e.vy),
         vx = e.vx + DT*f.fx,
         vy = e.vy + DT*f.fy,
         fx = f.fx,
         fy = f.fy
       }
 )
} ;;


// --- grow -------------------------------------------------------

DIFF  := 1.0  ;;
CONS  := 10.0  ;;
DiffP := 0.00023  ;; 
DeathSP := 0.00001 ;;
DivG := 0.00600 ;;
DEPOT := 15.0 ;; 
CroitG :=0.004  ;; 
CroitS :=0.007;;  


//cell division function
fun divide(b) = (
 b+{dap= b.dap/2, radius=R0_Gm },
 b + {dap= b.dap/2,x=noise(b.x),y=noise(b.y), radius=R0_Gm}
) ;;



trans Evol = {

 x / x.soma => if (random(1.0) <= DeathSP) then <undef> else if ( ( x.radius < R0_S ) & ( random(1.0) < CroitS) ) then x + 
{radius = x.radius + (R0_S - R0_G)/4} else  x fi fi ;  

 x => (
   let dap_diff = neighborsfold((\y.\acc.( {dap=DT*DIFF*(y.dap-x.dap) + acc.dap, n=acc.n+1} )), {dap=0.0,n=0}, x) in
   let dap' = if (x.dap + dap_diff.dap/dap_diff.n) - DT*CONS >= 0 then (x.dap + dap_diff.dap/dap_diff.n) - DT*CONS else 0 fi in

       if (random(1.0) <= DiffP)
       then x + {dap=DEPOT,soma=true,radius=x.radius} //la nouvelle cellule S a la taille de la cellule G dont elle provient
       else  if x.radius  >= R0_G then if  ( (random(1.0) <= DivG) & dap'> 1 )then divide(x + {dap=dap'}) else (x + {dap=dap'}) fi
 else  if random(1.0) < CroitG then x + {radius = x.radius + (R0_G - R0_Gm)/5 , dap=dap' } else x+{dap=dap'} fi fi fi
 ) 

} ;;




// --- Visualisation ----------------------------------------------------

outname := "/tmp/sheet.imo" ;;
outfile := (if (is_opened(outname)) then close(outname) else <undef> fi ; open(outname,1)) ;;




fun show_cell(e) = (
 "\tTranslated { Translation <" + e.x + ", " + e.y + ", " + 0.0 + "> Geometry Sphere { Radius " + e.radius + " Slices 16 
Stacks 16 " + "Color<" + if(e.soma) then 0 + ", " + 1.0 + ", " + 0 else if ((e.dap+0.4)*(0.5+R0_G/e.radius))< 1.0 then 
((e.dap+0.4)*(0.5+R0_G/e.radius)) else  1.0  fi+ ", " + 0 + ", " + 0 fi + "> } }"
) ;;



fun show[cpt=0](f, freq, c) = (
 cpt := 1 + cpt ;
 if (0 == cpt % freq)
 then (
   stdout << cpt << "\n" ;
   print_coll(f,
              c,
              show_cell,
              "Scene a"+cpt+" {\n",
              "\n",
              "}\nReplace { Show a"+cpt+"}\n\n")
 ) else <undef> fi ;
 c
) ;;

// --- Etat initial ----------------------------------------------------

v0 := {vx=0, vy=0, vz=0, fx=0, fy=0, fz=0, dap=14.0,soma=false,radius=R0_Gm} ;;
v1 :={vx=0, vy=0, vz=0, fx=0, fy=0, fz=0, dap=DEPOT,soma=true,radius=R0_S} ;;
pre_init :=
       v1+{x = 1.01,  y = 1.023, z = 0.101},
       v1+{x = 0.07,  y = 1.0,   z = 0.095},
       v1+{x = 1.0,   y = 0.01,  z = 0.098},
       v1+{x = 0.52,  y = 0.53,  z = 0.1},
       v1+{x = 0.01,  y = 0.02,  z = 0.099},
       v0+{x = -0.4,  y = 0.02,  z = 0.099},
       v0+{x = 1.4,  y = 0.02,  z = 0.099}
;;

init := Meca[iter=1000](delaunayfy(D2:(), pre_init)) ;;

//show(outfile, 1, init);;
//system("imoview "+outname);;

// --- Evolution -------------------------------------------------------

fun step(sys) = (
 Evol(Meca(sys))
) ;;

step[iter=100000,interlude=show(outfile,10)](init) ;;

system("./imoview_black "+outname);;
!quit ;;


Differentiation DAP dependent

record MecaBact = {x, y, vx, vy, fx, fy, radius}
and record CellBact  = {dap:float, soma:bool}
and record Bact = MecaBact + CellBact;;


delaunay(2) D2 = (\e.(e.x, e.y)) ;;

fun noise(x) = x + 0.0005 - random(0.001) ;;
// --- Mechanistic-------------------------------------------------------

DT    := 0.05;;
K     := 1.0;;
MU    := 1.8;;
R0_Gm := 0.50;; //germinal cells minmal rayon 
R0_G  :=0.75;;
R0_S  := 1.00;;

// interaction : the effect of the cells between each others
fun interaction(ref, src) = (
 let X      = ref.x - src.x
 and Y      = ref.y - src.y in
 let dist   = sqrt(X*X+Y*Y) in
 let spring = 0.0-K*(dist-(ref.radius+src.radius))/dist in
   {fx=X*spring - ref.vx*MU, fy = Y*spring - ref.vy*MU}
) ;;

fun add_vect(u, v) = { fx = u.fx + v.fx, fy = u.fy + v.fy } ;;
fun sum(x, u, acc) = add_vect(acc, interaction(x,u)) ;;


trans Meca = {
 e => (
   let f = neighborsfold(sum(e), {fx=0,fy=0}, e) in
     e+{ x  = noise(e.x  + DT*e.vx),
         y  = noise(e.y  + DT*e.vy),
         vx = e.vx + DT*f.fx,
         vy = e.vy + DT*f.fy,
         fx = f.fx,
         fy = f.fy
       }
 )
} ;;


// --- grow -------------------------------------------------------


DIFF  := 1.0  ;;
CONS  := 10.0  ;;
DiffP := 0.005  ;;
DeathSP := 0.00001 ;;
DivG := 0.00300 ;;
DEPOT := 16 ;;
CroitG :=0.002  ;; re
CroitS :=0.007;;  

fun divide(b) = (
 b+{dap= b.dap/2, radius=R0_Gm },
 b + {dap= b.dap/2,x=noise(b.x),y=noise(b.y), radius=R0_Gm}
) ;;



trans Evol = {

 x / x.soma => if (random(1.0) <= DeathSP) then <undef> else if ( ( x.radius < R0_S ) & ( random(1.0) < CroitS) ) then x + {radius
 = x.radius + (R0_S - R0_G)/4} else  x fi fi ; r)

  x => (
   let dap_diff = neighborsfold((\y.\acc.( {dap=DT*DIFF*(y.dap-x.dap) + acc.dap, n=acc.n+1} )), {dap=0.0,n=0}, x) in
   let dap' = if (x.dap + dap_diff.dap/dap_diff.n) - DT*CONS >= 0 then (x.dap + dap_diff.dap/dap_diff.n) - DT*CONS else 0 fi in
     if (dap'<=0.0)
     then (
       if (random(1.0) <= DiffP)
       then x + {dap=DEPOT,soma=true,radius=x.radius} //la nouvelle cellule S a la taille de la cellule G dont elle provient
       else x + {dap=dap'} fi
     ) else if x.radius  >= R0_G then if  ( (random(1.0) <= DivG) & dap'> 2 )then divide(x + {dap=dap'}) else (x + {dap=dap'}) fi 
else  if random(1.0) < CroitG then x + {radius = x.radius + (R0_G - R0_Gm)/5 , dap=dap' } else x+{dap=dap'} fi fi fi 
 ) 
} ;;




// --- Visualisation ----------------------------------------------------

outname := "/tmp/sheet.imo" ;;
outfile := (if (is_opened(outname)) then close(outname) else <undef> fi ; open(outname,1)) ;;




fun show_cell(e) = (
 "\tTranslated { Translation <" + e.x + ", " + e.y + ", " + 0.0 + "> Geometry Sphere { Radius " + e.radius + " Slices 16 
Stacks 16 " + "Color<" + if(e.soma) then 0 + ", " + 1.0 + ", " + 0 else if ((e.dap+0.4)*(0.5+R0_G/e.radius))< 1.0 then
 
((e.dap+0.4)*(0.5+R0_G/e.radius)) else  1.0  fi+ ", " + 0 + ", " + 0 fi + "> } }"
) ;;



fun show[cpt=0](f, freq, c) = (
 cpt := 1 + cpt ;
 if (0 == cpt % freq)
 then (
   stdout << cpt << "\n" ;
   print_coll(f,
              c,
              show_cell,
              "Scene a"+cpt+" {\n",
              "\n",
              "}\nReplace { Show a"+cpt+"}\n\n")
 ) else <undef> fi ;
 c
) ;;

// --- Etat initial ----------------------------------------------------

v0 := {vx=0, vy=0, vz=0, fx=0, fy=0, fz=0, dap=14.0,soma=false,radius=R0_Gm} ;;
v1 :={vx=0, vy=0, vz=0, fx=0, fy=0, fz=0, dap=DEPOT,soma=true,radius=R0_S} ;;
pre_init :=
       v1+{x = 1.01,  y = 1.023, z = 0.101},
       v1+{x = 0.07,  y = 1.0,   z = 0.095},
       v1+{x = 1.0,   y = 0.01,  z = 0.098},
       v1+{x = 0.52,  y = 0.53,  z = 0.1},
       v1+{x = 0.01,  y = 0.02,  z = 0.099},
       v0+{x = -0.4,  y = 0.02,  z = 0.099},
       v0+{x = 1.4,  y = 0.02,  z = 0.099}
;;

init := Meca[iter=1000](delaunayfy(D2:(), pre_init)) ;;

//show(outfile, 1, init);;
//system("imoview "+outname);;

// --- Evolution -------------------------------------------------------

fun step(sys) = (
 Evol(Meca(sys))
) ;;

step[iter=100000,interlude=show(outfile,10)](init) ;;

system("./imoview_black "+outname);;
!quit ;;

Gillespie Simulation

Global Variables

(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
(* 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.02 ;;  //0.2     ;; // Cre degradation
K_DAPiD  := 100000 ;;  //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.003 ;;  //0.05    ;; // Differentiation
K_DAPiP  := 0.6 ;;  //1.0     ;; // DAP production
K_DAPEx  := 1000000 ;; //0.5   ;; // DAP export
K_DAPIm  := 1e+09 ;;  //1.0     ;; // DAP import
K_DAPeD  := 20000 ;;  //0.1     ;; // DAP degradation

K_Div    := 0.0065 ;;  //0.02   ;;
K_Death  := 0.001 ;;  //0.01   ;;

nb_dape  := 0 ;;

fun globals_help() = (
  stdout << "\n** Help of the \"global\" module **\n\n" ;
  stdout << "Globals\n" ;
  stdout << "\t- Na (float): Avogadro number [6.0221415e23]\n\n" ;
  stdout << "\t- VolumeUnit (float): mean volume of a bacterium [1/Na]\n\n" ;
  stdout << "\t- K_CreD (float): Cre degradation reaction rate in bacteria [0.01]\n\n" ;
  stdout << "\t- K_DAPiD (float): DAP degradation reaction rate in bacteria [0.01]\n\n" ;
  stdout << "\t- K_DAPApI (float): DAPAp inhibition by DAP [2.0]\n\n" ;
  stdout << "\t- K_DAPApA (float): DAPAp re-activation [0.001]\n\n" ;
  stdout << "\t- K_CreP1 (float): Cre production with activated DAPAp [0.01]\n\n" ;
  stdout << "\t- K_CreP2 (float): Cre production with inhibited DAPAp [0.0]\n\n" ;
  stdout << "\t- K_Diff (float): bacteria differentition [0.1]\n\n" ;
  stdout << "\t- K_DAPiP (float): DAP production in bacteria [0.5]\n\n" ;
  stdout << "\t- K_DAPEx (float): DAP exportation [0.1]\n\n" ;
  stdout << "\t- K_DAPIm (float): DAP importation [0.1]\n\n" ;
  stdout << "\t- K_DAPeD (float): DAP degradation reaction rate in the environment [1.0]\n\n" ;

  stdout << "Functions\n" ;
  stdout << "\t- gamma(V) (float -> float):\n\t\tgamma function for reaction to stochastic constants translation\n\n"
) ;;

Structure

(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
(* State of the system                                                                            *)
(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
(* The system is composed of bacteria of different kinds: germinal and somatic                    *)
(*                                                                                                *)
(* It is modeled by nested multi-sets. The super multi-set corresponds to the environment where   *)
(* chemicals and bacteria are diffusing.                                                          *)
(* A set of functions is added to memorize the quantity of chemicals in the system in order not   *)
(* to evaluate it at each time.                                                                   *)
(* ********************************************************************************************** *)
(* ********************************************************************************************** *)

cptCre   := 0 ;;
cptDAPe  := 0 ;;
cptDAPi  := 0.0 ;;
cptDAP   := 0 ;;
cptBactS := 0 ;;
cptBactG := 0 ;;
cptBact  := 0 ;;
fun incrCre() = (cptCre := cptCre + 1) ;;
fun decrCre() = (cptCre := cptCre - 1) ;;
fun incrDAP() = (cptDAP := cptDAP + 1) ;;
fun decrDAP() = (cptDAP := cptDAP - 1) ;;
fun incrDAPe() = (cptDAPe := cptDAPe + 1) ;;
fun decrDAPe() = (cptDAPe := cptDAPe - 1) ;;
fun incrDAPi() = (cptDAPi := cptDAPi + 1) ;;
fun decrDAPi() = (cptDAPi := cptDAPi - 1) ;;
fun DAPe_DAPi() = (incrDAPi() ; decrDAPe()) ;;
fun DAPi_DAPe() = (incrDAPe() ; decrDAPi()) ;;
fun incrNCre(N) = (cptCre := cptCre + N) ;;
fun decrNCre(N) = (cptCre := cptCre - N) ;;
fun incrNDAP(N) = (cptDAP := cptDAP + N) ;;
fun decrNDAP(N) = (cptDAP := cptDAP - N) ;;
fun incrNDAPe(N) = (cptDAPe := cptDAPe + N) ;;
fun decrNDAPe(N) = (cptDAPe := cptDAPe - N) ;;
fun incrNDAPi(N) = (cptDAPi := cptDAPi + N) ;;
fun decrNDAPi(N) = (cptDAPi := cptDAPi - N) ;;

fun incrBact() = (cptBact := cptBact + 1) ;;
fun decrBact() = (cptBact := cptBact - 1) ;;
fun incrBactS() = (cptBactS := cptBactS + 1) ;;
fun decrBactS() = (cptBactS := cptBactS - 1) ;;
fun incrBactG() = (cptBactG := cptBactG + 1) ;;
fun decrBactG() = (cptBactG := cptBactG - 1) ;;
fun BactS_BactG() = (raise  (`Erreur "Impossible S to G differentiation.")) ;;
fun BactG_BactS() = (incrBactS() ; decrBactG()) ;;
fun incrNBact(N) = (cptBact := cptBact + N) ;;
fun decrNBact(N) = (cptBact := cptBact - N) ;;
fun incrNBactS(N) = (cptBactS := cptBactS + N) ;;
fun decrNBactS(N) = (cptBactS := cptBactS - N) ;;
fun incrNBactG(N) = (cptBactG := cptBactG + N) ;;
fun decrNBactG(N) = (cptBactG := cptBactG - N) ;;

collection Bact      = bag
and constraint BactG = [~`DAP_Box]Bact
and constraint BactS = [~`LOXP_Box && ~`LOXP_Box_Cre]Bact
and collection Env   = bag  ;;

id := 0 ;;
fun new_bactId() = (id := id+1 ; id) ;;
fun new_bact(t)  = (incrBact() ; Bactify({birth=t, date=0.0, id=new_bactId(), next_state=`Unknown}::seq:()) )  ;;
fun new_bactS(t) = (incrBactS() ; `DAPAp :: `DAP_Box  :: new_bact(t)) ;;
fun new_bactG(t) = (incrBactG() ; (*incrDAPi() ; incrDAP() ;*) `DAPAp :: `LOXP_Box :: new_bact(t)) ;;
fun new_env(t, nbG, nbS) = (
  incrNDAP(nb_dape) ;
  incrNDAPe(nb_dape) ;
  let env = fold((\i.\acc.((new_bactG(t))::acc)),Envify(fold((\n.\acc.(`DAP::acc)),seq:(),nb_dape)), nbG) in
    fold((\i.\acc.((new_bactS(t))::acc)),env, nbS)
) ;;
fun map_bact(f,env) = (
  map((\e.(if Bact(e) then f(e) else e fi)),env) //update_bact
) ;;

fun divide_bactG(t,b,id) = (
  let id' = new_bactId() in
  let nbDAP = count(`DAP,b)
  and nbCre = count(`Cre,b) in
  let bact1_core = Bactify({birth=t, date=0.0, id=id,  next_state=`Unknown}::seq:())
  and bact2_core = Bactify({birth=t, date=0.0, id=id', next_state=`Unknown}::seq:()) in
  let bact1_core = fold((\n.\acc.(`DAP::acc)),bact1_core,nbDAP/2)
  and bact2_core = fold((\n.\acc.(`DAP::acc)),bact2_core,nbDAP - nbDAP/2) in
  let bact2_core = fold((\n.\acc.(`Cre::acc)),bact2_core,nbCre/2)
  and bact1_core = fold((\n.\acc.(`Cre::acc)),bact1_core,nbCre - nbCre/2) in
  let bact1_core = if (member(`DAPAp_i,b)) then (`DAPAp_i) else (`DAPAp) fi :: bact1_core
  and bact2_core = `DAPAp :: bact2_core in
  let bact2_core = if (member(`LOXP_Box_Cre,b)) then (`LOXP_Box_Cre) else (`LOXP_Box) fi :: bact2_core
  and bact1_core = `LOXP_Box :: bact1_core in
    bact1_core :: bact2_core :: seq:()
) ;;

fun set_DAPe(N,env) = (
  let dDAP = cptDAPe - N in
    if (dDAP > 0)
    then (
      (* cptDAPe > N : trop de DAP *)
      cptDAPe := N ;
      cptDAP  := cptDAP - dDAP ;
      diff(env, fold((\n.\acc.(`DAP::acc)),Env:(),dDAP))
    )
    else (
      if (dDAP < 0)
      then (
	(* cptDAPe < N : pas assez de DAP *)
	cptDAPe := N ;
	cptDAP  := cptDAP - dDAP ;
	join(env, fold((\n.\acc.(`DAP::acc)),Env:(),-1*dDAP))
      )
      else (env) fi
    ) fi
) ;;









fun static_help() = (
  stdout << "\n** Help of the \"static\" module **\n\n" ;
  stdout << "Globals\n" ;
  stdout << "\t- cptCre (cptCre): `Cre counter\n\n" ;
  stdout << "\t- cptDAP (cptDAP): `DAP counter\n\n" ;
  stdout << "\t- cptDAPe (cptDAPe): `DAP counter (in the environment)\n\n" ;
  stdout << "\t- cptDAPi (cptDAPi): `DAP counter (in bacteria)\n\n" ;
  stdout << "\t- cptBact (cptBact): bacteria counter\n\n" ;
  stdout << "\t- cptBactG (cptBactG): germinal bacteria counter\n\n" ;
  stdout << "\t- cptBactS (cptBactS): somatic bacteria counter\n\n" ;

  stdout << "Functions\n" ;
  stdout << "\t- incrCre() (unit -> int):\n\t\tincrement cptCre\n\n" ;
  stdout << "\t- decrCre() (unit -> int):\n\t\tdecrement cptCre\n\n" ;
  stdout << "\t- incrDAP() (unit -> int):\n\t\tincrement cptDAP\n\n" ;
  stdout << "\t- decrDAP() (unit -> int):\n\t\tdecrement cptDAP\n\n" ;
  stdout << "\t- incrDAPe() (unit -> int):\n\t\tincrement cptDAPe\n\n" ;
  stdout << "\t- decrDAPe() (unit -> int):\n\t\tdecrement cptDAPe\n\n" ;
  stdout << "\t- incrDAPi() (unit -> int):\n\t\tincrement cptDAPi\n\n" ;
  stdout << "\t- decrDAPi() (unit -> int):\n\t\tdecrement cptDAPi\n\n" ;
  stdout << "\t- DAPe_DAPi() (unit -> int):\n\t\tdecrement cptDAPe, increment cptDAPi\n\n" ;
  stdout << "\t- DAPi_DAPe() (unit -> int):\n\t\tdecrement cptDAPi, increment cptDAPe\n\n" ;
  stdout << "\t- incrBact() (unit -> int):\n\t\tincrement cptBact\n\n" ;
  stdout << "\t- decrBact() (unit -> int):\n\t\tdecrement cptBact\n\n" ;
  stdout << "\t- incrBactG() (unit -> int):\n\t\tincrement cptBactG\n\n" ;
  stdout << "\t- decrBactG() (unit -> int):\n\t\tdecrement cptBactG\n\n" ;
  stdout << "\t- incrBactS() (unit -> int):\n\t\tincrement cptBactS\n\n" ;
  stdout << "\t- decrBactS() (unit -> int):\n\t\tdecrement cptBactS\n\n" ;
  stdout << "\t- BactG_BactS() (unit -> int):\n\t\tdecrement cptBactG, increment cptBactS\n\n" ;

  stdout << "\t- new_Bact(t) (float -> Bact):\n\t\tcreate a new bacterium at time t\n\n" ;
  stdout << "\t- new_BactG(t) (float -> BactG):\n\t\tcreate a new germinal bacterium at time t\n\n" ;
  stdout << "\t- new_BactS(t) (float -> BactS):\n\t\tcreate a new somatic bacterium at time t\n\n" ;
  stdout << "\t- new_env(t,ng,ns) (float * int * int -> Env):\n\t\tcreate a new environment containing ng germinal bacteria and ns somatic bacteria at time t\n\n" ;
  stdout << "\t- map_bact(f,env) ((Bact -> Bact) * Env -> Env):\n\t\tapply the function f on each bacterium of env\n\n" ;
  stdout << "\t- set_DAPe(N,env) (int * Env -> Env):\n\t\tchange env, cptDAP and cptDAPe to N `DAP\n\n"
) ;;

Scheduling

Dynamics

ExtSSA

Output for Gnuplot

Retrieved from "http://2007.igem.org/wiki/index.php/Paris/Sources"