Commit 1ddd2e36 authored by REMONDEAU Paul's avatar REMONDEAU Paul
Browse files

small changes to stoch conversion and comments

parent 82630201
......@@ -12,7 +12,8 @@
:- use_module(biocham).
:- use_module(numerical_simulation).
/* I don't know what to do with the stochastic conversion rate */
/* I don't know what to do with the stochastic conversion rate
It seems like species names like ABC-DFG-{HIJ} make the simulation not work, typical example MAPK*/
:- doc('This section describes a hybrid numerical simulation
method. Some reactions which occurs more often and which have
......@@ -155,6 +156,7 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
load_biocham(StochFilename),
all_items([kind: reaction],Reactions2),
nb_setval(reactions,Reactions2),
findall(R2,
(
......@@ -271,6 +273,8 @@ hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time)
write(StreamAppend,'option(show: {}).'),
close(StreamAppend),
nb_delete(volume),
nb_delete(reactions),
retractall(reaction_kinetics(_,_)),
load_biocham(OutFilename),
format(" ~n If you want to choose which species are plot, please change the corresponding ~n line in the output file (by default out_static.bc in the biocham directory). ~n ~n",[]).
......@@ -402,6 +406,7 @@ list_hybrid_ode(Stream,Hybrid_originals) :- % Only MA are converted, need to do
)
).
/* Can try to get the kinetics function in BIOCHAM 4 directly and translate, but i didn't know how
Also can direclty adapt the kinetics with the hybrid species list, but need to translate the kinetics first if needed */
......@@ -417,7 +422,7 @@ convert_MA(Reaction,Reactions,Hybrid_originals,NewKinetics) :-
atom_concat('k',ReactionNumber,Kdiff1),
atom_concat(Kdiff1,'_diff',Kdiff),
NewPartialKinetics = K * NewValue,
assertz(reaction(Reaction,NewPartialKinetics)),
assertz(reaction_kinetics(Reaction,NewPartialKinetics)),
NewKinetics = Kdiff * NewPartialKinetics
;
NewKinetics = K * NewValue
......@@ -473,7 +478,14 @@ convert_MA2(Reaction,ODESpecies,NewKinetics,Reactants,Products) :-
Kinetics =.. ['MA',K],
kinetics:eval_kinetics(Reactants, _, product('S'*'M' in [reactants], 'M'^'S'), Value),
convert_into_hybrid2(Value,ODESpecies,NewValue),!,
NewKinetics = K * NewValue.
NewKinetics = K * NewValue,
(
nb_current(dynamic_parameters,_)
->
true
;
assertz(reaction_kinetics(Reaction,NewKinetics))
).
%format("NewKinetics is ~w ~n",[NewKinetics]).
convert_into_hybrid2(Value * 1, ODESpecies,NewValue) :-
......@@ -499,9 +511,16 @@ convert_into_hybrid2(Molecule,ODESpecies,StochMolecule) :-
atom_concat(Molecule,'_stoch',StochMolecule),
!.
convert_into_hybrid2(Molecule,ODESpecies,HybridMolecule) :-
member(Molecule,ODESpecies),
atom_concat(Molecule,'_total',HybridMolecule),
convert_into_hybrid2(Molecule,Hybrid_originals,HybridMolecule) :-
get_option(stochastic_conversion, Rate),
nb_getval(volume,Volume),
member(Molecule,Hybrid_originals),
atom_concat('(',Molecule,Molecule1),
atom_concat(Molecule1,'_total/(',Molecule2),
atom_concat(Molecule2,Rate,Molecule3),
atom_concat(Molecule3,'*',Molecule4),
atom_concat(Molecule4,Volume,Molecule5),
atom_concat(Molecule5,'))',HybridMolecule),
!.
build_change_list(Changes, Counter, Length, Init, Result) :-
......@@ -549,7 +568,8 @@ summation(Input_list, Count, Init, Final) :-
)
).
print_event(Reactions_data, Stoch_species_order, _Stoch_reactants_order, Hybrid_reactants_list, Constraints_list, Stream) :-
print_event(Reactions_data, Stoch_species_order, _Stoch_reactants_order, Hybrid_reactants_list, Constraints_list, Stream) :-
get_option(stochastic_conversion,Rate),
length(Reactions_data, Number_of_alpha),
print_alpha_parameters(Number_of_alpha, 1, Stream),
print_constraints_parameters(Constraints_list, 1, Stream),
......@@ -557,7 +577,8 @@ print_event(Reactions_data, Stoch_species_order, _Stoch_reactants_order, Hybrid_
write(Stream, 'add_event(Time > tau,'),nl(Stream),
%%% print('Checkpoint 1\n'),
print_static_alpha(Number_of_alpha,Reactions_data, 1, Stream),
write(Stream,' tau = (if alpha_sum <= 0 then inf else Time+(-1 / alpha_sum) * log(random_float)),'),nl(Stream),
format("Checkpoint ~n",[]),
write(Stream,' tau = (if alpha_sum <= 0 then inf else Time+((-1 / alpha_sum) * log(random_float))/'),write(Stream,Rate),write(Stream,'),'),nl(Stream),
write(Stream,' ran = random_float,'),nl(Stream),
%print_event_constraints(Constraints_list, 1, Stream),
%concat_species(Stoch_species_order, 1, '', Species_for_event),
......@@ -650,10 +671,11 @@ print_static_alpha(Number_of_alpha,Reactions_data, Current_count, Stream) :-
write(Stream, ' alpha'),
write(Stream, Current_count),
write(Stream, ' = '),
nth1(Current_count,Reactions_data,Reaction),
nth1(2,Reaction,Kinetics),
convert_kinetics(Kinetics,NewKinetics1),
simplify(NewKinetics1,NewKinetics),
nb_getval(reactions,Reactions),
nth1(Current_count,Reactions,Reaction),
reaction_kinetics(Reaction,NewKinetics),
format("Reaction ~w is ~w Kinetics is ~w ",[Current_count,Reaction,NewKinetics]),
(
nb_current(dynamic_parameters,_)
->
......@@ -1288,6 +1310,7 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,ParticuleTreshold,PropTresh, SortedStochSpecies, StochReactants,Hybrid_reactants_list, Constraints_list, StreamAppend),
nb_delete(dynamic_parameters),
nb_delete(reactions),
retractall(reaction_kinetics(_,_)),
nl(StreamAppend),
write(StreamAppend,'% Plot part'),nl(StreamAppend),nl(StreamAppend),
......@@ -1304,6 +1327,11 @@ hybrid_dynamic_simulation(InputFile,OutFilename,Volume,Time,PropTresh,PartTresh)
load_biocham(OutFilename),
format(" ~n If you want to choose which species are plot, please change the corresponding ~n line in the output file (by default out_dynamic.bc in the biocham directory). ~n ~n",[]).
:- doc('\\begin{example}T7 model using hybrid dynamic simulation:').
:- biocham_silent(clear_model).
:- biocham(hybrid_dynamic_simulation('library:examples/hybrid/Goutsias_dynamic.bc',20,10)).
:- doc('\\end{example}').
compute_partial_maximum_change(RSum,PSum) :-
retract(minimum_reactant_change(OldRMin)),
retract(maximum_reactant_change(OldRMax)),
......@@ -1405,6 +1433,7 @@ print_dynamic_parameters(NumberOfReactions,Stream) :-
!.
print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,MaxParticleChange,PropTresh,Stoch_species_order, _Stoch_reactants_order, Hybrid_reactants_list, Constraints_list, Stream) :-
get_option(stochastic_conversion,Rate),
length(Reactions_data, Number_of_alpha),
print_alpha_parameters(Number_of_alpha, 1, Stream),
print_constraints_parameters(Constraints_list, 1, Stream),
......@@ -1415,7 +1444,7 @@ print_dynamic_event(SortedSpecies,Reactions_data, HybridDynamicReactants,MaxPart
print_dynamic_event_condition(SortedSpecies,HybridDynamicReactants,1,PropTresh,MaxParticleChange,Stream),
print_static_alpha(Number_of_alpha,Reactions_data, 1, Stream),
write(Stream,' delta_t = (if alpha_sum <= 0 then step_size else (-1 / alpha_sum) * log(random_float)),'),nl(Stream),
write(Stream,' tau = Time+delta_t,'),nl(Stream),
write(Stream,' tau = Time+delta_t/'),write(Stream,Rate),write(Stream,','),nl(Stream),
write(Stream,' ran = random_float,'),nl(Stream),
%print_event_constraints(Constraints_list, 1, Stream),
%concat_species(Stoch_species_order, 1, '', Species_for_event),
......
MA(c1) for gen => tem.
MA(c2) for tem => _.
MA(c3) for tem => tem + gen.
MA(c4) for gen + struc => virus.
MA(c5) for tem => tem + struc.
MA(c6) for struc => _.
parameter(c1 = 0.025).
parameter(c2 = 0.25).
parameter(c3 = 1).
parameter(c4 = 0.0000075).
parameter(c5 = 1000).
parameter(c6 = 1.99).
present(gen,0).
present(tem).
present(struc,0).
present(virus,0).
\ No newline at end of file
MA(k1) for A + 2*B => C.
MA(k2) for C => 2*A.
parameter(k1 = 1).
parameter(k2 = 1).
present(A,1).
present(B,1).
present(C,0).
\ No newline at end of file
parameter(tau = 0).
parameter(ran = 0).
parameter(delta_t = 0).
parameter(k1_diff = 0).
parameter(k2_diff = 0).
function(kk1 = k1_diff).
function(kk2 = k2_diff).
function(A_total = floor(100*1*A + A_stoch)).
parameter(A_stoch = 0).
function(AA = A_stoch).
function(B_total = floor(100*1*B + B_stoch)).
parameter(B_stoch = 0).
function(BB = B_stoch).
function(C_total = floor(100*1*C + C_stoch)).
parameter(C_stoch = 0).
function(CC = C_stoch).
k1_diff*(k1*((A_total/(100*1))*(B_total/(100*1))^2))for A+2*B=>C.
k2_diff*(k2*(C_total/(100*1)))for C=>2*A.
present(A,1).
present(B,1).
present(C,0).
parameter(
k1 = 1,
k2 = 1
).
parameter(alpha1 = 0).
parameter(alpha2 = 0).
parameter(alpha_sum = 0).
parameter(lower_level1 = 0).
parameter(lower_level2 = 0).
parameter(lower_level3 = 0).
add_event(Time > tau,
k1_diff = (if k1*(A_total*B_total^2) >= 1/(step_size * 20) and A_total >= 10 and B_total >= 10 then 1 else 0),
k2_diff = (if k2*C_total >= 1/(step_size * 20) and C_total >= 20 then 1 else 0),
alpha1 = (1-k1_diff)*(k1/2*A_total/100*B_total/100*(B_total+ -1)/100),
alpha2 = (1-k2_diff)*(k2*C_total/100),
alpha_sum = alpha1 + alpha2,
delta_t = (if alpha_sum <= 0 then step_size else (-1 / alpha_sum) * log(random_float)),
tau = Time+delta_t/100,
ran = random_float,
lower_level1 = (if A_total >= 1 then 1 else 0),
lower_level2 = (if B_total >= 2 then 1 else 0),
lower_level3 = (if C_total >= 1 then 1 else 0),
A_stoch = (if alpha_sum*ran>0 and alpha_sum*ran<=alpha1 and lower_level1>0 then A_stoch-1
else if alpha_sum*ran>alpha1 and alpha_sum*ran<=alpha1 + alpha2 and lower_level2>0 then A_stoch+2
else A_stoch),
B_stoch = (if alpha_sum*ran>0 and alpha_sum*ran<=alpha1 and lower_level1>0 then B_stoch-2
else B_stoch),
C_stoch = (if alpha_sum*ran>0 and alpha_sum*ran<=alpha1 and lower_level1>0 then C_stoch+1
else if alpha_sum*ran>alpha1 and alpha_sum*ran<=alpha1 + alpha2 and lower_level2>0 then C_stoch-1
else C_stoch)
).
% Plot part
option(show: {A_total, B_total, C_total}). % Choose here the species to plot
numerical_simulation(method:rsbk, time: 20).
plot.
option(show: {}).
\ No newline at end of file
MA(r2a__a2) for P_KKK+E2=>E2_P_KKK.
MA(r2a__d2) for E2_P_KKK=>P_KKK+E2.
MA(r2b__k2) for E2_P_KKK=>E2+KKK.
MA(r3a__a3) for KK+P_KKK=>P_KKK_KK.
MA(r3a__d3) for P_KKK_KK=>KK+P_KKK.
MA(r3b__k3) for P_KKK_KK=>P_KK+P_KKK.
MA(r4a__a4) for P_KK+KKPase=>KKPase_P_KK.
MA(r4a__d4) for KKPase_P_KK=>P_KK+KKPase.
MA(r4b__k4) for KKPase_P_KK=>KK+KKPase.
MA(r5a__a5) for P_KK+P_KKK=>P_KKK_P_KK.
MA(r5a__d5) for P_KKK_P_KK=>P_KK+P_KKK.
MA(r5b__k5) for P_KKK_P_KK=>PP_KK+P_KKK.
MA(r6a__a6) for PP_KK+KKPase=>KKPase_PP_KK.
MA(r6a__d6) for KKPase_PP_KK=>PP_KK+KKPase.
MA(r6b__k6) for KKPase_PP_KK=>P_KK+KKPase.
MA(r7a__a7) for K+PP_KK=>PP_KK_K.
MA(r7a__d7) for PP_KK_K=>K+PP_KK.
MA(r7b__k7) for PP_KK_K=>P_K+PP_KK.
MA(r8a__a8) for P_K+KPase=>KPase_P_K.
MA(r8a__d8) for KPase_P_K=>P_K+KPase.
MA(r8b__k8) for KPase_P_K=>K+KPase.
MA(r9a__a9) for P_K+PP_KK=>PP_KK_P_K.
MA(r9a__d9) for PP_KK_P_K=>P_K+PP_KK.
MA(r9b__k9) for PP_KK_P_K=>PP_KK+PP_K.
MA(r10a__a10) for PP_K+KPase=>KPase_PP_K.
MA(r10a__d10) for KPase_PP_K=>PP_K+KPase.
MA(r10b__k10) for KPase_PP_K=>P_K+KPase.
present(E2,0.0003).
present(KKK,0.003).
present(P_KKK,0.0).
present(KK,1.2).
present(P_KK,0.0).
present(PP_KK,0.0).
present(K,1.2).
present(P_K,0.0).
present(PP_K,0.0).
present(KPase,0.12).
present(KKPase,0.0003).
present(E2_P_KKK,0.0).
present(P_KKK_KK,0.0).
present(P_KKK_P_KK,0.0).
present(PP_KK_K,0.0).
present(PP_KK_P_K,0.0).
present(KKPase_PP_KK,0.0).
present(KKPase_P_KK,0.0).
present(KPase_PP_K,0.0).
present(KPase_P_K,0.0).
parameter(
compartment = 4.0e-12,
K_PP_norm_max = 0.900049,
r2a__a2 = 1000.0,
r2a__d2 = 150.0,
r2b__k2 = 150.0,
r3a__a3 = 1000.0,
r3a__d3 = 150.0,
r3b__k3 = 150.0,
r4a__a4 = 1000.0,
r4a__d4 = 150.0,
r4b__k4 = 150.0,
r5a__a5 = 1000.0,
r5a__d5 = 150.0,
r5b__k5 = 150.0,
r6a__a6 = 1000.0,
r6a__d6 = 150.0,
r6b__k6 = 150.0,
r7a__a7 = 1000.0,
r7a__d7 = 150.0,
r7b__k7 = 150.0,
r8a__a8 = 1000.0,
r8a__d8 = 150.0,
r8b__k8 = 150.0,
r9a__a9 = 1000.0,
r9a__d9 = 150.0,
r9b__k9 = 150.0,
r10a__a10 = 1000.0,
r10a__d10 = 150.0,
r10b__k10 = 150.0
).
function(
K_PP_norm = (PP_K+KPase_PP_K)/(PP_K+P_K+K+PP_KK_K+PP_KK_P_K+KPase_PP_K+KPase_P_K),
rel_K_PP_max = K_PP_norm/K_PP_norm_max,
KK_PP_norm = (PP_KK+PP_KK_K+PP_KK_P_K+KKPase_PP_KK)/(PP_KK+P_KK+KK+PP_KK_K+PP_KK_P_K+P_KKK_KK+P_KKK_P_KK+KKPase_PP_KK+KKPase_P_KK),
KKK_P_norm = (P_KKK+P_KKK_KK+P_KKK_P_KK)/(KKK+P_KKK+P_KKK_KK+P_KKK_P_KK)
).
\ No newline at end of file
MA(r1a__a1) for KKK+E1=>E1_KKK.
MA(r1a__d1) for E1_KKK=>KKK+E1.
MA(r1b__k2) for E1_KKK=>E1+P_KKK.
present(E1,3.0e-5).
present(KKK,0.003).
present(E1_KKK,0.0).
present(P_KKK,0.0).
parameter(
r1a__a1 = 1000.0,
r1a__d1 = 150.0,
r1b__k2 = 150.0).
\ No newline at end of file
......@@ -26,8 +26,12 @@ stochastic_simulation(Until, Propensities, Time, Boolean) :-
retractall(last_time(_)),
assertz(last_time(0.0)),
get_option(stochastic_thresholding, Threshold),
%statistics(runtime,_),
markov_run(InitState, Conditions, Actions, Until, Propensities, Time,
Threshold),
%statistics(runtime,[_,T]),
%TT is T/1000,
%format("Simulation time: ~ps~n",[TT]),
make_state_table,
restore_parameters.
......
......@@ -377,7 +377,7 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
nb_setval(k_time_units,2)
),
repeat,
%counter(TTTTT), write(TTTTT), nl,
counter(TTTTT), write(TTTTT), nl,
nb_getval(parameters_list,Parameters),
nb_getval(numerical_table,CurrentTable),
get_last_row(CurrentTable,Time,CurrentStateInitial),
......@@ -910,7 +910,7 @@ scalar_multiply([Element|Vector],Scalar,[NewElement|Result]) :-
scalar_multiply(Vector,Scalar,Result).
euler(Equations,Time,CurrentState,CurrentTable) :-
%format("One event occured ~n",[]),
format("One event occured ~n",[]),
retractall(euler_last_state(_,_)),
assertz(euler_last_state(Time,CurrentState)),
nb_getval(euler_step_size,EulerStepSize),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment