Paris/Sources
From 2007.igem.org
(Difference between revisions)
(→Constant rate of differentiation) |
(→Differentiation DAP dependent) |
||
| Line 274: | Line 274: | ||
| - | // --- | + | // --- Growth ----------------------------------------------------- |
Revision as of 15:32, 25 October 2007
Contents |
Cellular Automaton for DAP Diffusion
!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) ;;
Cells Spatial Organization
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
}
)
} ;;
// --- Growth ------------------------------------------------------
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
}
)
} ;;
// --- Growth -----------------------------------------------------
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 Based 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.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.5 ;; // DAP export
K_DAPIm := 1.0 ;; // DAP import
K_DAPeD := 0.1 ;; // DAP degradation
K_Div := 0.02 ;; // Mitosis rate
K_Death := 0.01 ;; // Death rate
nb_dape := 0 ;; // Number of DAP molecule in the environment (initial condition)
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
constraint schedule = `Empty_SC | schedule_entry
and record schedule_entry = { date:float, id:int, sc:schedule, rule:int } ;;
fun schedule_empty() = `Empty_SC ;;
fun schedule_add_entry(id,date,rule,sc) = (
switch (sc)
case `Empty_SC : {date=date,id=id,rule=rule,sc=sc}
default : (
if (sc.date > date)
then {date=date,id=id,rule=rule,sc=sc}
else {date=sc.date,id=sc.id,rule=sc.rule,sc=schedule_add_entry(id,date,rule,sc.sc)} fi
)
endswitch
) ;;
fun schedule_remove_entry(id,sc) = (
switch (sc)
case `Empty_SC : `Empty_SC
default : (
if (sc.id == id)
then sc.sc
else {date=sc.date,id=sc.id,rule=sc.rule,sc=schedule_remove_entry(id,sc.sc)} fi
)
endswitch
) ;;
fun schedule_size(sc) = (
switch (sc)
case `Empty_SC : 0
default : 1 + schedule_size(sc.sc)
endswitch
) ;;
fun schedule_replace_entry(id,date,rule,sc) = (
switch (sc)
case `Empty_SC : {date=date,id=id,rule=rule,sc=sc}
default : (
if (sc.id == id)
then schedule_add_entry(id,date,rule,sc.sc)
else (
if (sc.date > date)
then {date=date,id=id,rule=rule,sc=schedule_remove_entry(id,sc)}
else {date=sc.date,id=sc.id,rule=sc.rule,sc=schedule_replace_entry(id,date,rule,sc.sc)} fi
) fi
)
endswitch
) ;;
fun schedule_remove_first_entry(sc) = sc.sc ;;
fun schedule_print(sc) = (
switch (sc)
case `Empty_SC : stdout << "\n"
default : (
stdout << sc.id << "\t(" << sc.rule << ", " << sc.date << ") \n" ;
schedule_print(sc.sc)
)
endswitch
) ;;
schedule := schedule_empty() ;;
fun schedule_help() = (
stdout << "\n** Help of the \"schedule\" module **\n\n" ;
stdout << "Globals\n" ;
stdout << "\t- schedule (schedule): the main schedule initilized with value `Empty_SC\n\n" ;
stdout << "Functions\n" ;
stdout << "\t- schedule_empty() (unit -> schedule):\n\t\treturns the empty schedule\n\n" ;
stdout << "\t- schedule_add_entry(id,d,r,sc) (int * float * int * schedule -> schedule):\n\t\tadd an entry in the schedule sc corresponding to the\n\t\tapplication of the rule r on the bacterium id at the date d\n\n" ;
stdout << "\t- schedule_add_entry(id,d,r,sc) (int * float * int * schedule -> schedule):\n\t\treplace (or add if not existing) an entry in the\n\t\tschedule sc corresponding to the application of the rule r on the bacterium id at the date d\n\n" ;
stdout << "\t- schedule_remove_entry(id,sc) (int * schedule -> schedule):\n\t\tremove in sc the first occurence of an event in bacterium id\n\n" ;
stdout << "\t- schedule_remove_first(sc) (schedule -> schedule):\n\t\tremove the first entry of sc\n\n" ;
stdout << "\t- schedule_print(sc) (schedule -> stdout):\n\t\tprint sc\n\n"
) ;;
Dynamics
(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
(* Dynamics of the system *)
(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
(* It exists two kinds of dynamics. The first is related to the chemical reactions that occur *)
(* everywhere in the system (inside and outside the bacteria + transport of proteins between the *)
(* outside and inside). This parts will be simulated using a sotchastic process based on the *)
(* Gillespie's algorithm. The second part of the dynamics to be modeled corresponds to the *)
(* bacteria behavior like divisions and cells death. In fact, a lot of processes appearing in *)
(* bacteria and belonging to their normal dynamics are skipped/abstracted as functions of the *)
(* chemical state of bacteria. For example, a too high concentration of DAP is lethal ; we *)
(* reprensent this dependance by a function that deletes all bacteria where the threshold is *)
(* exceeded. These behavior functions have to be evaluated after each chemical update. *)
(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
constraint ruleId = `CreD | `DAPiD | `DAPApI | `DAPApA | `CreP1 | `CreP2 | `Diff1 | `Diff2 | `DAPiP | `DAPEx | `DAPIm | `DAPeD | `Div | `Death | `Unknown ;;
fun fstOrder(k,v) = k ;;
fun sndOrder(k,v) = k/gamma(v) ;;
//fun count(p,c) = fold((\e.\acc.(acc + if p(e) then 1 else 0 fi)),0,c) ;;
//fun countAll(p1,p2,c) = fold((\e.\acc.(acc + if p2(e) then count(p1,e) else 0 fi)),0,c) ;;
appliedBactRule := `Unkown ;;
trans inBact = {
(* 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 = \c.(fstOrder(K_CreP1* if (member(`DAP,c)) then 1.0 else 1.0 fi,1.0)) }=> (appliedBactRule := `CreP1 ; `DAPAp, `Cre) ;
`DAPAp_i ={ C = fstOrder(K_CreP2,1.0) }=> (appliedBactRule := `CreP2 ; `DAPAp_i, `Cre) ;
(* Differentiation *)
`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 *)
`DAP ={ C = fstOrder(K_DAPEx,1.0) }=> (appliedBactRule := `DAPEx ; return(`DAP :: (diff(self,`DAP)) :: seq:())) ;
} ;;
appliedEnvRule := `Unkown ;;
nbReleasedDAP := 0 ;;
deadBactId := -1 ;;
trans inEnv[upBact,current_t,current_sc] = {
(* DAP internalization *)
`DAP, b:Bact ={ A = (\c.(sndOrder(K_DAPIm,1.0) * cptDAPe * cptBact)) }=> (appliedEnvRule := `DAPIm ; (upBact(current_sc,(`DAP :: b))) ) ;
(* DAP external degradation *)
`DAP ={ C = fstOrder(K_DAPeD,1.0) }=> (appliedEnvRule := `DAPeD ; <undef>) ;
(* BactS death *)
(x:BactS \/ {id=id}) as b ={ A = (\c.(fstOrder(K_Death,1.0)*cptBactS)) }=> (
let nbDAP = count(`DAP,b) in
nbReleasedDAP := nbDAP ;
appliedEnvRule := `Death ;
deadBactId := id ;
map((\e.(`DAP)),iota(nbDAP,seq:()))
) ;
(* BactG mitosis *)
(_:BactG \/ ({id=id})) as b ={ A = (\c.(fstOrder(K_Div,1.0)*cptBactG)) }=> (
appliedEnvRule := `Div ;
let b1b2 = divide_bactG(current_t+'tau,b,id) in
let b1 = upBact(current_sc,b1b2.(0)) in
let b2 = upBact(schedule,b1b2.(1)) in
b1 :: (b2 :: seq:())
) ;
} ;;
fun cptUpdate[nbRD](r) = ( // Stochiometric update for the output
switch (r)
case `CreD: decrCre()
case `DAPiD: (decrDAPi() ; decrDAP())
case `DAPApI: (decrDAPi() ; decrDAP())
case `DAPApA: (incrDAPi() ; incrDAP())
case `CreP1: incrCre()
case `CreP2: incrCre()
case `Diff1: decrCre()
case `Diff2: (decrCre(); BactG_BactS())
case `DAPiP: (incrDAPi() ; incrDAP())
case `DAPEx: DAPi_DAPe()
case `DAPIm: DAPe_DAPi()
case `DAPeD: (decrDAPe() ; decrDAP())
case `Death: (decrNDAPi(nbRD) ; incrNDAPe(nbRD) ; decrBactS() ; decrBact())
case `Div: (incrBactG() ; incrBact())
default: raise (`Erreur ("Diet_coli: Not yet implemented: rule "+r))
endswitch
) ;;
fun dynamic_help() = (
stdout << "\n** Help of the \"dynamic\" module **\n\n" ;
stdout << "Requirements\n" ;
stdout << "\t- module \"globals\"\n" ;
stdout << "\t- module \"static\"\n\n" ;
stdout << "Globals\n" ;
stdout << "\t- appliedBactRule (ruleId): takes the id of a rule after the application of inBact [`Ùnknown]\n\n" ;
stdout << "\t- appliedEnvRule (ruleId): takes the id of a rule after the application of inEnv [`Ùnknown]\n\n" ;
stdout << "Functions\n" ;
stdout << "\t- fstOrder(k,v) (float * float -> float):\n\t\tkinetics to stochastic conversion\n\n" ;
stdout << "\t- sndOrder(k,v) (float * float -> float):\n\t\tkinetics to stochastic conversion\n\n" ;
stdout << "\t- sndOrder(k,v) (float * float -> float):\n\t\tkinetics to stochastic conversion\n\n" ;
stdout << "\t- inBact(b) (Bact -> Bact | `DAP*Bact):\n\t\tchemical reactions in bacteria ; to be applied with strategy `gillespie\n\n" ;
stdout << "\t- inEnv[update_bact(Bact -> Bact)](e) (Env -> Env):\n\t\tchemical reactions in environment ; the option update_bact is to update the state of a bacterium modified by the application of the rule\n\n" ;
stdout << "\t- cptUpdate(r) (ruleId -> unit):\n\t\tupdate of the global counters from \"static\" w.r.t. the applied rule r\n\n"
) ;;
ExtSSA
(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
(* System evolution algorithm *)
(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
(* The model is a one depth level P system. Our alogrithm consists in applying the gillespie SSA *)
(* on each bacterium. A bacterium contains the current state together with the next step and the *)
(* corresponding reaction delay. The SSA is applied on the environment providing a reaction time. *)
(* The evolution is done where the delay is the smallest. *)
(* ********************************************************************************************** *)
(* ********************************************************************************************** *)
trans no_next_state = { { next_state } as x => x + { next_state = `Unknown, date=0.0 } } ;;
trans update_bact[t,sc] = {
{ id=id } as x => (
let date = 0.0 in
let nx_state = inBact[strategy=`gillespie,
postlude=(\c.(date := t+'tau ; c)),
iter=1](no_next_state(self))
in
schedule := schedule_replace_entry(id,date,appliedBactRule,sc) ;
appliedBactRule := `Unkown ;
x + { date=date, next_state = nx_state }
)
} ;;
trans evolve_bact[t,sc,idEv] = {
((b:Bact \/ ({id=x,next_state=nxs}))) as bb / (idEv == x) => (
if (Bact(nxs))
then update_bact[t=t,sc=sc](nxs)
else (nxs.(0)::(update_bact[t=t,sc=sc](nxs.(1)) :: seq:())) fi
)
} ;;
fun step[get_global_time,set_global_time,sc](env) = (
t := get_global_time() ;
let date' = 0.0 in
let env' = inEnv[upBact=(\sc.\x.(let t' = t+'tau in update_bact[t=t',sc=sc,strategy=`default](x))),
current_t = t,
current_sc = sc,
strategy=`gillespie,
postlude=(\c.(date' := t + 'tau ; c)),
iter=1](env)
in
//!! (stdout << "schedule.date=" << schedule.date << " > date'=" << date' << "\n" ; (date' == infinity) || (date' < schedule.date)) ;
if (date' < sc.date)
then (
//stdout << "envionment evolution (applied rule = " << appliedEnvRule << "(" << nbReleasedDAP << "), date = " << date' << "(" << t << "))\n" ;
cptUpdate[nbRD=nbReleasedDAP](appliedEnvRule) ;
schedule := if (appliedEnvRule==`Death) then schedule_remove_entry(deadBactId,schedule) else schedule fi ;
nbReleasedDAP := 0 ;
deadBactId := -1 ;
appliedEnvRule := `Unkown
)
else (
//stdout << "bacterium evolution (id = " << sc.id << ", applied rule = " << sc.rule << ", date = " << sc.date << "(" << t << "))\n" ;
cptUpdate(sc.rule) ;
nbReleasedDAP := 0 ;
deadBactId := -1 ;
appliedEnvRule := `Unkown ;
date' := sc.date ;
env' := evolve_bact[strategy=`asynchronous,sc=sc,idEv=sc.id,t=date'](env)
) fi ;
//stdout << env' << "\n" ;
sc := schedule ;
set_global_time(date') ;
env'
) ;;
fun evol_algo_help() = (
stdout << "\n** Help of the \"evol_algo\" module **\n\n" ;
stdout << "Requirements\n" ;
stdout << "\t- module \"schedule\"\n" ;
stdout << "\t- module \"dynamic\"\n\n" ;
stdout << "Functions\n" ;
stdout << "\t- update_bact[t=t(float),sc=sc(schedule)](b) (Bact -> Bact):\n\t\tcompute the next Gillespie step of a bacterium b. The schedule sc is updated and the new value of sc is returned in the global variable schedule. The option t has to be set to the global date.\n\n" ;
stdout << "\t- evolve_bact[t=t(float),sc=sc(schedule),idEv=id(int)](env) (Env -> Env):\n\t\tmake the bacterium with identifiant id evolve. The schedule sc is updated and the new value of sc is returned in the global variable schedule. The option t has to be set to the global date.\n\n" ;
stdout << "\t- step[t=t(float),sc=sc(schedule)](env) (Env -> Env):\n\t\tOne step evolution of the whole system. The earliest event is applied in the environment env w.r.t. the schedule sc. The option t has to be set to the global date.\n\n"
) ;;
Output for Gnuplot
then (
let tg = get_global_time() in
GPoutput << get_global_time() << "\t" <<
cptCre << "\t" <<
cptDAP << "\t" <<
cptDAPe << "\t" <<
(cptDAPi/cptBact) << "\t" <<
cptBact << "\t" <<
cptBactS << "\t" <<
cptBactG << "\n" ;
stdout << "Iteration: " << 'iteration << " and Global time: " << tg << " and bacteria population: " << cptBact << "\n"
)
else <undef> fi ;
env
) ;;
fun dataToPng() = (
gnuplot << "set key outside\n";
gnuplot << "set xlabel \"Time (arbitrary unit)\"\n";
gnuplot << "set ylabel \"Number of molecules or structures\"\n";
gnuplot << "set output \""+ PNGoutput + "\"\n";
gnuplot << "set term png\n";
gnuplot << "plot " +
"\""+ GPoutput +"\" using 1:4 title 'DAPe' with lines 3," +
"\""+ GPoutput +"\" using 1:5 title 'DAPi' with lines 4," +
"\""+ GPoutput +"\" using 1:6 title 'Bact' with lines 5," +
"\""+ GPoutput +"\" using 1:7 title 'BactS' with lines 6," +
"\""+ GPoutput +"\" using 1:8 title 'BactG' with lines 7" +
"" ;
gnuplot << "\nquit" ;
gnuplot << "\n"
) ;;
fun dataToPngDAP() = (
gnuplot << "set key outside\n";
gnuplot << "set xlabel \"Time (arbitrary unit)\"\n";
gnuplot << "set ylabel \"Number of molecules or structures\"\n";
gnuplot << "set output \""+ PNGoutput + "\"\n";
gnuplot << "set term png\n";
gnuplot << "plot " +
"\""+ GPoutput +"\" using 1:4 title 'DAPe' with lines 3," +
"\""+ GPoutput +"\" using 1:5 title 'DAPi' with lines 4" +
"" ;
gnuplot << "\nquit" ;
gnuplot << "\n"
) ;;
fun dataToPngBact() = (
gnuplot << "set key outside\n";
gnuplot << "set xlabel \"Time (arbitrary unit)\"\n";
gnuplot << "set ylabel \"Number of molecules or structures\"\n";
gnuplot << "set output \""+ PNGoutput + "\"\n";
gnuplot << "set term png\n";
gnuplot << "plot " +
"\""+ GPoutput +"\" using 1:6 title 'Bact' with lines 5," +
"\""+ GPoutput +"\" using 1:7 title 'BactS' with lines 6," +
"\""+ GPoutput +"\" using 1:8 title 'BactG' with lines 7" +
"" ;
gnuplot << "\nquit" ;
gnuplot << "\n"
) ;;
fun output_help() = (
stdout << "\n** Help of the \"output\" module **\n\n" ;
stdout << "Requirements\n" ;
stdout << "\t- module \"static\"\n" ;
stdout << "Globals\n" ;
stdout << "\t- GPoutput (string): Gnuplot formatted output file [\"/tmp/sortie.dat\"]\n\n" ;
stdout << "\t- PNGoutput (string): PNG file for Gnuplot exporting (see function dataToPng) [\"/tmp/sortie.png\"]\n\n" ;
stdout << "\t- gnuplot (string): absolute path to Gnuplot executable [\"/usr/bin/gnuplot\"]\n\n" ;
stdout << "Functions\n" ;
stdout << "\t- export_cpt[get_global_time](freq, env) (int * Env -> Env):\n\t\tExport in GPoutput the state of global module counters (to be used as interlude of an iterated function)\n\n" ;
stdout << "\t- dataToPng() (unit -> unit):\n\t\texport GPoutput to the PNG file PNGoutput using Gnuplot (located at path gnuplot)\n\n"
) ;;
Main program
!set match_random := 100 ;;
duration := 1500.0 ;;
gt := 0.0 ;; //global time
fun ggt() = gt ;; //global time set and get function... an assert is used to check that time only increases
fun sgt(x) = (!! (gt <= x) ; gt := x) ;;
!include "globals.mgs" ;;
!include "static.mgs" ;;
!include "schedule.mgs" ;;
!include "dynamic.mgs" ;;
!include "evol_algo.mgs" ;;
!include "output.mgs" ;;
schedule := schedule_empty() ;;
env := map_bact((\b.(update_bact[t=gt,sc=schedule](b))),new_env(gt,10,0)) ;;
step[fixpoint=(\c1.\c2.(gt > duration || cptBact == 0)), interlude=B(export_cpt[get_global_time=ggt](10),I(*set_DAPe(nb_dape)*)), sc=schedule, get_global_time=ggt, set_global_time=sgt](env) ;;
dataToPng() ;;
