Paris/Sources
From 2007.igem.org
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) ;;
Gillespie Simulation
Global Variables
Structure
Scheduling
Dynamics
ExtSSA
Output for Gnuplot
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 ;;
