Commit 0537fa7f authored by FAGES Francois's avatar FAGES Francois
Browse files

notebooks ASSB finaux

parent 7ba0e0d0
......@@ -363,10 +363,10 @@ satisfaction_degree(Formula, Objective) :-
type(Objective, [variable_objective]),
option(foltl_magnitude, number, _Magnitude, 'order of magnitude for << and >> operators'),
doc('
computes a continuous satisfaction degree in the interval [0,+∞) for a
FOLTL(Rlin) constraint on the current trace. The degree is greater or
equal than 1 if the formula is satisfied. The greater the degree the more
robust the satisfaction of the formula.
computes a continuous satisfaction degree in the interval [0,+∞) for an
FOLTL(Rlin) property on the current trace with respect to some objective values for the free variables of the formula. The degree is greater or
equal than 1 if the formula is satisfied. The greater the degree the greater the margin in the satisfaction of the formula (i.e. formula satisfaction robustness).
This satisfaction degree is computed by (an approximation of) the distance of the objective point to the validity domain of the formula (if less than 1), or by its penetration depth (if greater than 1).
'),
satisfaction_degree(Formula, Objective, Degree),
format('~f\n', [Degree]).
......
%% Cell type:markdown id: tags:
# ADVANCES IN SYSTEMS AND SYNTHETIC BIOLOGY
## Modelling Complex Biological Systems in the Context of Genomics
## Thematic Research School
https://assb.lri.fr 20 March 2018
# Workshop on the BIOCHAM modeling environment
François Fages and Sylvain Soliman, Inria Saclay Ile de France,
[ASSB BIOCHAM workshop](https://assb.lri.fr/en/Abstracts/BIOCHAM.html)
%% Cell type:markdown id: tags:
# Minimal mitotic oscillator
after *[Albert Goldbeter 1991 PNAS](https://doi.org/10.1073/pnas.88.20.9107)*
<img src="oscillator.png" width=300>
%% Cell type:code id: tags:
```
load(oscillator.bc).
```
%% Cell type:code id: tags:
```
list_model.
```
%% Cell type:code id: tags:
```
draw_reactions.
```
%% Cell type:markdown id: tags:
**Question:** _Do you think this system can oscillate? Does it have a_ **negative feedback loop**?
%% Cell type:code id: tags:
```
draw_influences.
```
%% Cell type:markdown id: tags:
**Question:** _Can you find which species are (mass-)conserved in this system?_
%% Cell type:code id: tags:
```
search_conservations.
```
%% Cell type:markdown id: tags:
## Differential semantics
Ordinary differential equations (ODEs) on continuous molecular concentrations
%% Cell type:code id: tags:
```
list_ode.
```
%% Cell type:code id: tags:
```
numerical_simulation. plot.
```
%% Cell type:code id: tags:
```
plot(show:Cyclin, against:ProteaseP).
```
%% Cell type:markdown id: tags:
**Question:** _Can you find parameter values that increase or decrease Cyclin's period?_
%% Cell type:code id: tags:
```
%slider vi kd V2 K2
```
%% Cell type:markdown id: tags:
## Stochastic semantics
Continuous time Markov Chain (CTMC) on integer numbers of molecules (conversion factor 100 for concentration 1)
%% Cell type:code id: tags:
```
parameter(vi=0.025, kd=0.01, V2=1.5, K2=0.005).
```
%% Cell type:code id: tags:
```
numerical_simulation(method:ssa). plot.
```
%% Cell type:markdown id: tags:
## Boolean semantics
Non-deterministic asynchronous Boolean transition system on Boolean states of present/absent molecules
## Computation Tree Logic (CTL) for querying the Boolean transition paths.
Logical connectives: /\ (and), \/ (or), not (negation), -> (implication)
Path quantifiers: E (there exists a path), A (for all paths)
Time operators: F (finally at some time point), G (globally at all time points), U (until)
Example: reachable state EF(s), steady state EG(s), stable state AG(s), etc.
%% Cell type:code id: tags:
```
generate_ctl.
```
%% Cell type:code id: tags:
```
expand_ctl(reachable(s)).
```
%% Cell type:code id: tags:
```
expand_ctl(steady(s)).
```
%% Cell type:code id: tags:
```
expand_ctl(stable(s)).
```
%% Cell type:code id: tags:
```
expand_ctl(oscil(f)).
```
%% Cell type:markdown id: tags:
## Model reduction preserving CTL properties
Which reactions can be removed while preserving the (generated) CTL specification of the behavior ?
%% Cell type:code id: tags:
```
reduce_model.
```
%% Cell type:code id: tags:
```
numerical_simulation. plot.
```
%% Cell type:code id: tags:
```
load(oscillator.bc).
```
%% Cell type:code id: tags:
```
numerical_simulation. plot.
```
%% Cell type:markdown id: tags:
## Linear Time Logic with Constraints over Reals FO-LTL(Rlin)
No path quantifier (single trace analysis)
Only temporal operators F, G, U, X (next)
Free variables over the reals
Linear constraints over real variables
Existential (exists) and universal (forall) quantifiers over real variables
%% Cell type:code id: tags:
```
validity_domain(amplitude(Cyclin,v)).
```
%% Cell type:code id: tags:
```
expand_ltl(amplitude(x,a)).
```
%% Cell type:code id: tags:
```
validity_domain(period(Cyclin,p)).
```
%% Cell type:code id: tags:
```
satisfaction_degree(period(Cyclin,p), [p->25]).
```
%% Cell type:code id: tags:
```
satisfaction_degree(period(Cyclin,p), [p->22]).
```
%% Cell type:markdown id: tags:
## Parameter search
%% Cell type:code id: tags:
```
seed(0). search_parameters(period(Cyclin,p), [0<=vi<=1, 0<=V2<=2], [p -> 22]).
```
%% Cell type:code id: tags:
```
list_parameters.
```
%% Cell type:code id: tags:
```
numerical_simulation. plot.
```
%% Cell type:code id: tags:
```
validity_domain(period(Cyclin,p)).
```
%% Cell type:markdown id: tags:
---
Let's now have a closer look at the amplitude of the _Cyclin_…
%% Cell type:code id: tags:
```
plot(show: Cyclin).
```
%% Cell type:code id: tags:
```
validity_domain(amplitude(Cyclin, amp)).
```
%% Cell type:markdown id: tags:
**Question:** _Does this formula capture the amplitude between peaks?_
_Using the_ `local_minimum(Species, min)` _and_ `local_maximum(Species, max)` _constructs, can you write a pattern to capture that?_
%% Cell type:code id: tags:
```
expand_ltl(local_minimum(x,m)).
```
%% Cell type:code id: tags:
```
ltl_pattern(peak_amplitude(Species, Amplitude) =
exists(min, exists(max, local_minimum(Species, min) /\ local_maximum(Species, max) /\ Amplitude <= max - min))).
```
%% Cell type:code id: tags:
```
validity_domain(peak_amplitude(Cyclin, amp)).
```
%% Cell type:code id: tags:
```
seed(0). search_parameters(peak_amplitude(Cyclin, amp), [0<=vi<=1, 0<=V2<=2], [amp -> 0.8]).
```
%% Cell type:code id: tags:
```
option(show: {Cyclin, Protease, Kinase}).
```
%% Cell type:code id: tags:
```
numerical_simulation. plot.
```
%% Cell type:code id: tags:
```
validity_domain(peak_amplitude(Cyclin, amp)).
```
%% Cell type:markdown id: tags:
## Robustness of the model w.r.t. an FO-LTL(Rlin) property and parameter perturbations
* Model robustness (mean satisfaction degree for normal distribution of the parameter values)
* Formula satisfaction robustness optimization (maximization of the penetration depth in the validity domain)
%% Cell type:code id: tags:
```
robustness(peak_amplitude(Cyclin, amp), [vi, V2], [amp -> 0.8]).
```
%% Cell type:code id: tags:
```
seed(0). search_parameters(peak_amplitude(Cyclin, amp), [0<=vi<=1, 0<=V2<=2], [amp -> 0.8], cmaes_stop_fitness: -0.1).
```
%% Cell type:code id: tags:
```
numerical_simulation. plot.
```
%% Cell type:code id: tags:
```
validity_domain(peak_amplitude(Cyclin, amp)).
```
%% Cell type:code id: tags:
```
robustness(peak_amplitude(Cyclin, amp), [vi, V2], [amp -> 0.8]).
```
%% Cell type:code id: tags:
```
validity_domain(period(Cyclin, p)).
```
%% Cell type:markdown id: tags:
## Synthesis of continuous Biochemical Reaction Networks for Real Functions
* Definition of computable real functions by polynomial initial value problems (PIVPs)
* Synthesis of continuous CRNs from PIVPs (showing Turing completeness of PIVPs)
Synthesis of CRNs for the cosine function as a function of time or of an input x
%% Cell type:code id: tags:
```
clear_model.
```
%% Cell type:code id: tags:
```
compile_from_expression(cos,time,f).
```
%% Cell type:code id: tags:
```
list_model.
```
%% Cell type:code id: tags:
```
option(time:10). numerical_simulation. plot.
```
%% Cell type:code id: tags:
```
numerical_simulation(method:ssa). plot.
```
%% Cell type:code id: tags:
```
clear_model.
```
%% Cell type:code id: tags:
```
compile_from_expression(cos,x,f).
```
%% Cell type:code id: tags:
```
list_model.
```
%% Cell type:code id: tags:
```
present(x_p,4).
```
%% Cell type:code id: tags:
```
numerical_simulation. plot.
```
%% Cell type:code id: tags:
```
seed(0). numerical_simulation(method:ssa). plot.
```
%% Cell type:code id: tags:
```
dose_response(x_p, 0, 6.28, show:{f_p, f_m}).
```
%% Cell type:markdown id: tags:
## Subsidiary question
Implement the cos(time) function with activation reactions (e.g. by phosphorylation) instead of synthesis reactions
Propose a deactivation scheme for the annihilation reactions.
%% Cell type:code id: tags:
```
```
......
This diff is collapsed.
......@@ -12,6 +12,90 @@
:- use_module(objects).
:- use_module(biocham).
:- doc('The continuous satisfaction degree of an FO-LTL(Rlin) in a given trace with respect to the objective values for the free variables can be used to compute parameter sensitivity indices and robustness measures with respect to parameter perturbation according to normal distributions, and to search parameter values for satisfying an FO-LTL(Rlin) formula or even maximizing the margins.').
:- initial('option(robustness_samples: 100)').
:- initial('option(robustness_relative_error: 0.05)').
:- initial('option(robustness_coeff_var: 0.1)').
robustness(Formula, Parameters, Objective) :-
biocham_command,
type(Formula, foltl),
type(Parameters, [parameter_name]),
type(Objective, [variable_objective]),
option(time, time, _Time, 'time horizon of the numerical integration'),
option(method, method, _Method, 'method for the numerical solver'),
option(
robustness_coeff_var,
number,
Variation,
'coefficient of variation of the normal law'
),
option(
robustness_samples,
integer,
Samples,
'number of samples used for averaging'
),
option(
robustness_relative_error,
number,
RelativeError,
'relative sampling error used to stop estimation of the robustness'
),
option(
error_epsilon_absolute, number, _ErrorEpsilonAbsolute,
'absolute error for the numerical solver'
),
option(
error_epsilon_relative, number, _ErrorEpsilonRelative,
'relative error for the numerical solver'
),
option(
initial_step_size, number, _InitialStepSize,
'initial step size for the numerical solver'
),
option(
maximum_step_size, number, _MaximumStepSize,
'maximum step size for the numerical solver, as a fraction of time'
),
option(
precision, number, _Precision,
'precision for the numerical solver'
),
doc('
computes the robustness degree as defined in \\cite{RBFS11tcs}, with
respect to formula \\argument{Formula}, for the list of parameters
\\argument{Parameters} and with list of objectives for the free variables
of \\argument{Formula} given in \\argument{Objective}.
The robustness is the (sampled) average from a normal law with
coefficient of variation (stddev/mean) given by the corresponding option
and centered around the current parameter values where at each point the
satisfaction degree is computed
(and truncated to 1 if necessary).
\\clearmodel
\\begin{example}
\\trace{
biocham: k for a => b. present(a,10). parameter(k=0.7).
biocham: robustness(G(Time < T => a > b), [k], [T -> 5]).
}
\\end{example}
'),
with_current_ode_system(
with_clean(
[
numerical_simulation:variable/2,
numerical_simulation:equation/2,
numerical_simulation:parameter_index/2,
numerical_simulation:conditional_event/2
],
robustness_aux(Parameters, [(Formula, [], Objective)], Samples,
RelativeError, Variation)
)
).
:- initial('option(cmaes_init_center: no)').
:- initial('option(cmaes_log_normal: no)').
......@@ -57,19 +141,16 @@ search_parameters(Formula, Parameters, Objective) :-
),
option(
cmaes_stop_fitness, number, _Fitness,
'stop when distance to the objective is less than this value (use
negative to enforce robustness)'
'stop when distance to the objective is less than this value (use 0 by default and
negative values in [-1,0] to maximize margins)'
),
doc('
tries to satisfy a FOLTL(Rlin) constraint by varying the parameters listed
tries to satisfy a FO-LTL(Rlin) constraint by varying the parameters listed
in \\argument{Parameters}.
This search command uses the stochastic optimization procedure CMA-ES
(covariance matrix adaptation evolution strategy)
with the continuous satisfaction degree of a FO-LTL(Rlin) property as
with the continuous satisfaction degree (truncated to 1 or not according to stop option) of the given FO-LTL(Rlin) property as
fitness function.
This satisfaction degree is defined by the distance between the validity
domain of a FOLTL(Rlin) formula with free variables (as given by the command
domains) and the objective values given in \\argument{Objective}.
\\clearmodel
\\begin{example}
\\trace{
......@@ -138,8 +219,7 @@ search_parameters(Parameters, Conditions) :-
negative to enforce robustness)'
),
doc('similar to \\command{search_parameters/3} but uses the list of
\\argument{Conditions} as triples with an FOLTL formula, a mutant
described as a list of parameters instantiations, and an objective for the
\\argument{Conditions} as triples with an FOLTL formula, a list of parameters instantiations (e.g. describing a mutant), and an objective for the
variables of the formula.
\\clearmodel
\\begin{example}'),
......@@ -165,88 +245,6 @@ search_parameters(Parameters, Conditions) :-
list_ids(ParameterIds).
:- initial('option(robustness_samples: 100)').
:- initial('option(robustness_relative_error: 0.05)').
:- initial('option(robustness_coeff_var: 0.1)').
robustness(Formula, Parameters, Objective) :-
biocham_command,
type(Formula, foltl),
type(Parameters, [parameter_name]),
type(Objective, [variable_objective]),
option(time, time, _Time, 'time horizon of the numerical integration'),
option(method, method, _Method, 'method for the numerical solver'),
option(
robustness_coeff_var,
number,
Variation,
'coefficient of variation of the normal law'
),
option(
robustness_samples,
integer,
Samples,
'number of samples used for averaging'
),
option(
robustness_relative_error,
number,
RelativeError,
'relative sampling error used to stop estimation of the robustness'
),
option(
error_epsilon_absolute, number, _ErrorEpsilonAbsolute,
'absolute error for the numerical solver'
),
option(
error_epsilon_relative, number, _ErrorEpsilonRelative,
'relative error for the numerical solver'
),
option(
initial_step_size, number, _InitialStepSize,
'initial step size for the numerical solver'
),
option(
maximum_step_size, number, _MaximumStepSize,
'maximum step size for the numerical solver, as a fraction of time'
),
option(
precision, number, _Precision,
'precision for the numerical solver'
),
doc('
computes the robustness degree as defined in \\cite{RBFS11tcs}, with
respect to formula \\argument{Formula}, for the list of parameters
\\argument{Parameters} and with list of objectives for the free variables
of \\argument{Formula} given in \\argument{Objective}.
The robustness is the (sampled) average from a normal law with
coefficient of variation (stddev/mean) given by the corresponding option
and centered around the current parameter values where at each point the