Chapter 1
Overview
1.1. Main features
The Biochemical Abstract Machine (Biocham) is a modelling software for cell systems biology, with some unique features for static and dynamic analyses using temporal logic constraints.
Biocham is a free software protected by the GNU General Public License GPL
version 2. This is an Open Source license that allows free usage of this
software.
This reference manual (as its extended version for developers) is automatically generated from the source code of Biocham.
Biocham v4 is mainly composed of :

a rulebased language for modeling biochemical reaction networks and influence networks
 compatible with SBML2 [doi.org/10.1038/npre.2008.2715.1] and SBML3qual [doi.org/10.1186/175205097135] respectively
 interpreted in a hierarchy of semantics (differential, stochastic, Petri Net, Boolean)
 a temporal logic based language to formalize possibly imprecise behaviours both qualitative and quantitative [RBFS11tcs]}
 several unique features for developing⁄analyzing⁄correcting⁄completing⁄calibrating⁄coupling⁄synthesizing reaction and influence networks with respect to formal specifications of their behaviour.
Biocham v4 is a complete rewriting (in SWIProlog) of Biocham v3 (written in GNUProlog).
It introduces some new features, mainly:
 influence network models [FMRS18tcbb]
 quantitative temporal logic patterns and trace simplifications [FS14book]
 compilation of mathematical functions and simple programs in reaction networks [FLBP17cmsb]
 PAC learning of influence models from data [CFS17cmsb]
 a notebook based on Jupyter
plus since v4.1:
 multistationarity check [BFS18jtb]}
 robustness optimization [FS18cmsb]
 detection of model reductions by constrained subgraph epimorphisms (SEPI) [GSF10bi]
 tropical equilibrations [SFR14amb]
 a menu of commands in the notebook
 a Graphical User Interface based on the notebook
Some features of Biocham v3 are still not implemented in v4, mainly:
 the graphical editor (SBGN compatible)
 hybrid stochasticdifferential simulations [CFJS15tomacs]
or less efficiently (because of onthefly C compilation instead of native Prolog code):
 numerical integration using GSL library is slower and may be less accurate
 parameter search, sensitivity and robustness measures are currently slower than Biocham v3.
1.2. Biocham files, notebooks and call options
Biocham files
Biocham file names are suffixed by
.bc
. In a Biocham file, everything following the character percent is a comment.
Some other file formats are used.
Biocham models can be imported from, and exported to, other file formats using the following suffixes:

.xml
for Systems Biology Markup Language (SBML) files, more precisely SBML2 files for reaction networks and SBML3qual files for influence networks; 
.ode
for Ordinary Differential Equation files (ODEs in XXPAUT format), allowing us to infer and import a reaction network from ODEs using a heuristic inference algorithm described in [FGS15tcs]. 
.ipynb
for Jupyter notebook files.
Biocham numerical data time series can be imported⁄exported in

.csv
commaseparated values format (spreadsheet format).
The following files can also be used to export some Biocham objects:

.tex
LaTeX format for exporting ODEs and graphs; 
.dot
for exporting graphs; 
.plot
for exporting numerical data time series; 
.smv
for exporting boolean transition systems and Computation Tree Logic (CTL) queries for the NuSMV modelchecker; 
.dot
for exporting graphs.
Biocham call options
biocham [file]
launches Biocham and optionnally loads the file given as argument.biocham –trace –version –listcommands
are possible options (prefixed by two dashes).biocham_debug
is the command to use to launch Biocham with the Prolog toplevel.Biocham kernel Jupyter notebook
biocham –notebook [notebook.ipynb]
launches Jupyter (see
http://jupyternotebook.readthedocs.io/) with a biocham kernel, and
loads the corresponding notebook (ipynb suffix).
If no notebook is provided,
a biocham notebook can be created with the Jupyter menu "new".
The
–lab
option will try to launch JupyterLab instead of the
traditional notebook.–console
will use the terminal console
of Jupyter.
The environment variable
BIOCHAM_TIMEOUT
can be used to change the
default notebookkernel communication timeout of 120 (seconds).
To execute a Biocham command enter ShiftReturn.
All shortcuts are described in the keyboard menu.
The following video
https://www.youtube.com/watch?v=jZ952vChhuI is a quite nice
introduction to the Jupyter notebook's concepts, all can be adapted to
Biocham as a kernel.
The only magic commands available in the notebook and not in Biocham are
%load file.bc
which creates a cell for each command contained in a Biocham file,
and%slider k1 k2 …
, which creates sliders to change the given
parameters' value, and run numerical_simulation
followed by
plot
at each change.
Biocham GUI
The graphical user interface of Biocham4 is based on the Jupyter notebook.
It provides a static view of the current state with buttons for all Biocham commands.
A button allows you to switch between the notebook timeline view and the GUI spatial view.
1.3. Interpreter toplevel
The exhaustive list of Biocham commands available at toplevel is described in the rest of this manual.
New reactions and new influences can also be entered directly at toplevel, as a short hand for the commands
add_reaction/1
and add_influence/1
respectively.
Some commands (e.g.,
numerical_simulation/0
in the example below) take named options
as arguments.
All the options can be defined either locally for a single command or
globally for the whole model.
Options can be defined for a single command
by adding additional arguments of the form Option: Value
.
Options can be defined globally for the model with the command
option(Option: Value, ..., Option: Value)
.
Local options take precedence over global options.
list_options.
lists the set of options defined in the current model with their default values.
Example.
biocham:
list_options.From inherited 'initial': [0] option(draw_first:present) [1] option(left_to_right:yes) [2] option(force_graph:no) [3] option(output_to_library:no) [4] option(import_reactions_with_inhibitors:no) [5] option(conservation_size:4) [6] option(mapping_restriction:[]) [7] option(merge_restriction:no) [8] option(timeout:180) [9] option(all_reductions:no) [10] option(distinct_species:no) [11] option(max_nb_reductions:200) [12] option(extremal_sepi:no) [13] option(michaelis_menten:yes) [14] option(r_1:yes) [15] option(r_2:no) [16] option(ep:no) [17] option(enzyme:yes) [18] option(hill_reaction:no) [19] option(partial_hill_reaction:no) [20] option(double_michaelis_menten:no) [21] option(michaelis_menten_expansion:no) [22] option(tropical_epsilon:0.1) [23] option(tropical_max_degree:3) [24] option(tropical_ignore:[{}]) [25] option(tropical_denominator:0) [26] option(tropical_single_solution:no) [27] option(test_permutations:no) [28] option(test_transpositions:no) [29] option(boolean_semantics:negative) [30] option(nusmv_initial_states:all) [31] option(nusmv_counter_example:no) [32] option(query:current_spec) [33] option(data_transform:none) [34] option(cnf_clause_size:3) [35] option(boolean_simulation:no) [36] option(method:bsimp) [37] option(error_epsilon_absolute:1.0e6) [38] option(error_epsilon_relative:1.0e6) [39] option(initial_step_size:1.0e6) [40] option(maximum_step_size:0.01) [41] option(minimum_step_size:1.0e5) [42] option(precision:6) [43] option(time:20) [44] option(stochastic_conversion:100) [45] option(stochastic_thresholding:1000) [46] option(filter:no_filter) [47] option(stats:no) [48] option(logscale:'') [49] option(against:Time) [50] option(xmin:auto) [51] option(ymin:auto) [52] option(xmax:auto) [53] option(ymax:auto) [54] option(foltl_magnitude:5) [55] option(robustness_samples:100) [56] option(robustness_relative_error:0.05) [57] option(robustness_coeff_var:0.1) [58] option(cmaes_init_center:no) [59] option(cmaes_log_normal:no) [60] option(cmaes_stop_fitness:0.0001) [61] option(simplify_variable_init:no) [62] option(zero_ultra_mm:no) [63] option(fast_rate:1000) [64] option(binomial_reduction:no) [65] option(lazy_negatives:yes) [66] option(show:[{}])
biocham:
option(time:100).Example.
biocham:
a=>b.biocham:
b=>a.biocham:
present(a).biocham:
list_ode.[0] d(a)/dt=ba [1] d(b)/dt=ab
biocham:
option(time:5).biocham:
numerical_simulation(method:ssa). plot.option.
sets options globally.
quit.
quits the interpreter.
seed(number).
For the sake of reproducibility of the results with randomized algorithms, sets the seed of the random number generator to some integer. Note however than the results may still differ on different machines.
prolog(Term: name).
Just for development purposes, calls a Prolog term written between quotes and prints the result substitution.
keep_temporary_files.
Switch to a mode where Biocham commands keep the temporary files they
generate.
delete_temporary_files.
Switch to a mode where Biocham commands delete the temporary files they
generate. (Default)
reset_options.
Reset all options to their default values.
about.
lists the version, copyright, license and URL of Biocham
Part I: Biochemical Networks
This part describes the syntax of BIOCHAM models of biochemical processes. They are essentially of two forms:
 reaction networks (i.e. Feinberg's Chemical Reaction Networks, CRNs, compatible with SBML)
 influence networks (variant of Thomas's regulatory networks, compatible with qualSBML).
Both types of models can be combined and interpreted in a hierarchy of semantics, including the differential semantics (Ordinary Differential Equations), stochastic semantics (Continuoustime Markov Chain), Petri net semantics, and Boolean semantics including several variants defined by options [FMRS18tcbb].
Chapter 2
Biochemical Objects
2.1. Syntax
The syntax of Biocham is described with formal grammar rules which define new syntactic tokens from primitive tokens such as atom (i.e. string), number, term (e.g. atom(..., ...)).
For instance, the syntax of an input or output file is just the syntax of an atom in both cases, but they are distinguished in this manual for documentation purposes:
input_file ::= 
output_file ::= 
Note that these atoms can enclose with single quotes (') usual UNIX
file expansion wildcards like * or ?.
The syntax of elementary types in Biocham is the following:
concentration ::= 
time ::= 
name ::= 
parameter_name ::= 
function_prototype ::= 
object ::= 
arithmetic_expression ::= 
condition ::= 
 true  false 
yesno ::= 
 yes  no 
2.2. Molecules and initial state
list_molecules.
lists all the molecules of the current model.
When a molecule is written as compound@location, it represents the
given compound as located inside the compartment location
list_locations.
lists all the locations of the current model.
list_neighborhood.
lists all neighborhood relationships inferred from the model.
draw_neighborhood.
draws the graph of all neighborhood relationships inferred from the
model.
list_initial_state.
lists the objects which are present (including their initial concentration)
and absent from the initial state.
clear_initial_state.
makes undefined all objects possibly present or absent in the initial
state.
present({object_{1}, ...,
object_{n}}).
Every object in
[object_{1}, ...,
object_{n}]
is made present in the initial
state with initial concentration equal to 1.
present({object_{1}, ...,
object_{n}}, concentration).
Every object in
[object_{1}, ...,
object_{n}]
is initialized
with the given initial concentration
, which can be a parameter
name to indicate that the parameter value will be used.
An initial value equal to 0 is equivalent to absent.
absent({object_{1}, ...,
object_{n}}).
undefined({object_{1}, ...,
object_{n}}).
Every object in
[object_{1}, ...,
object_{n}]
is made possibly present or absent
in the initial state. This is the default.
make_present_not_absent.
makes all objects (appearing in the instances of the current set of rules)
which are not declared absent, present in the initial state.
make_absent_not_present.
makes all objects (appearing in the instances of the current set of rules)
which are not declared present, absent in the initial state.
2.3. Parameters
parameter(parameter_name_{1} = number_{1},
...,
parameter_name_{n} =
number_{n}).
sets the value of parameters.
list_parameters.
shows the values of all known parameters.
delete_parameter(parameter_name_{1},
...,
parameter_name_{n}).
deletes some parameters
2.4. Functions
function(function_prototype_{1} = arithmetic_expression_{1},
...,
function_prototype_{n} =
arithmetic_expression_{n}).
sets the definition of functions.
list_functions.
lists all known functions.
delete_function(functor_{1},
...,
functor_{n}).
deletes some functions. Either arity is given, or all functions with
the given functor are deleted.
2.5. Aliases
alias(object_{1} = ... =
object_{n}).
makes
Objects
be alternative names for the
same object.
biocham:
a=>b.biocham:
b=>c.biocham:
c=>a.biocham:
list_reactions.[0] MA(1) for a=>b [1] MA(1) for b=>c [2] MA(1) for c=>a
biocham:
alias(a=c).biocham:
list_reactions.[0] MA(1) for a=>b [1] MA(1) for b=>a [2] MA(1) for a=>a
canonical(object).
makes
object
be the canonical representant
for all its aliases.
list_aliases.
shows the values of all known aliases.
delete_alias(object).
makes
object
distinct from all other objects.
Chapter 3
Reaction Networks
(Bio)chemical reaction networks (CRNs) can be described in BIOCHAM by a reaction model,
i.e. a finite (multi)set of directed reaction rules, possibly given with a rate function, using the syntax defined by the grammar below.
BIOCHAM reaction models are compatible with the Systems Biology Markup Language SBML version 2.
reaction ::= 
rule_name ::= 
basic_reaction ::= 
reactants ::= 
kinetics ::= 
Useful abbreviations for mass action law kinetics (with inhibitors), MichaelisMenten kinetics, Hill kinetics (with inhibitors).
function MA(k) = k * product(S * M in [reactants], M ^ S).
function MAI(k) = k * product(S * M in [reactants], M ^ S) / (1 + sum(M in [inhibitors], M)).
function MM(Vm, Km) = Vm * single_reactant / (Km + single_reactant)
.
function Hill(Vm, Km, n) =
Vm * single_reactant ^ n / (Km ^ n + single_reactant ^ n).
.
function HillI(Vm, Km, n) = Vm * single_reactant ^ n / (Km ^ n + single_reactant ^ n + sum(M in [inhibitors], M ^ n))..
3.1. Reaction editor
add_reaction(reaction).
adds reaction rules to the current set of reactions.
This command is implicit: reaction rules are directly added when reading them.
Remark. In Biocham v4, a reaction network is represented by a multiset of reactions. Reactions can thus be added in multiple copies in which case their reaction rates are summed.
delete_reaction(reaction).
removes one reaction rule from the current set of reactions.
delete_reaction_named(name).
removes reaction rules by their name from the current set of reactions.
delete_reactions.
Removes all reaction rules from the current model.
delete_reaction_prefixed(Prefix: name).
removes the reaction rules that match a name prefix.
find_reaction_prefixed(Prefix: name, name).
merge_reactions(reaction_{1}, reaction_{2}).
merges two reaction rules by removing them and replacing them by a new reaction with reactants the sum of reactants, with products the sum of the products and with the sum of the kinetics.
list_reactions.
lists the current set of reaction rules.
list_rules.
lists the current set of reaction and influence rules.
3.2. Reaction graph
up_type ::= 
 none  present  sources 
Defines if in the drawing of a graph using Graphviz, some species should appear first (i.e. on top of the drawing, or at the left if dran from left to right. The default value is
present
, i.e. draw first the species that are present in the initial state.option(draw_first:present).
The bipartite graph of molecular species and reactions can be drawn with the following command.
draw_reactions.
Draws the reaction graph of the current model.
Equivalent to
reaction_graph. draw_graph.
Example.
biocham:
load(library:examples/mapk/mapk).biocham:
draw_reactions.reaction_graph.
Builds the reaction graph of the current model.
Example.
biocham:
option(draw_first:none,left_to_right:no).biocham:
draw_reactions.Options.
 draw_first: up_type
 Put species of this type at the top of the graph when drawing it.
rule_graph.
Builds the rule graph of the current model, i.e. the union of
the reaction graph and the influence graph.
import_reactions_from_graph.
Updates the set of reactions of the current model with the current graph.
draw_rules.
Draws the reaction graph of the current model.
Equivalent to
rule_graph. draw_graph.
3.3. Graph visualization and editing
edgeref ::= 
new_graph.
Creates a new graph.
delete_graph(name).
Deletes a graph.
set_graph_name(name).
Sets the name of the current graph.
list_graphs.
Lists the graph of the current model.
select_graph(name).
Selects a graph
add_vertex(name_{1},
...,
name_{n}).
Adds new vertices to the current graph.
delete_vertex(name_{1},
...,
name_{n}).
Deletes a set of vertices from the current graph.
add_edge(edge_{1},
...,
edge_{n}).
Adds the given set of edges to the current graph.
The vertices are added if needed.
delete_edge(edgeref_{1},
...,
edgeref_{n}).
Deletes a set of edges from the current graph.
list_edges.
Lists the edges of the current graph.
list_isolated_vertices.
Lists the isolated vertices of the current graph.
list_graph_objects.
Lists the edges and the isolated vertices of the current graph,
and their attributes if any.
set_attribute({graph_object_{1}, ...,
graph_object_{n}}, attribute).
Adds an attribute to every vertex or edge in the given set.
The vertices and the edges are added if needed.
place(name_{1},
...,
name_{n}).
transition(name_{1},
...,
name_{n}).
delete_attribute(graph_object, Attribute: name).
Removes an attribute from
graph_object
.
list_attributes(graph_object).
List the attributes of
graph_object
.
draw_graph.
Draws the current graph.
Options.
 left_to_right: yesno
 Draws the graph from left to right instead of the default top to bottom.
option(left_to_right:yes).
export_graph(output_file).
Exports the current graph in a file.
The format is chosen from the suffix:
.dot
, .pdf
, .eps
, .ps
,
.png
or .svg
– assuming no extension is .dot
.
Chapter 4
Influence Networks
BIOCHAM models can also be defined by influence systems with forces, possibly mixed to reactions with rates.
In the syntax described by the grammar below, one influence rule (either positive ">" or negative "<") expresses that a conjunction of sources (or their negation after the separator "⁄" for inhibitors) has an influence on a target molecular species. This syntax necessitates to write the Boolean activation (resp. deactivation) functions of the molecular species in Disjunctive Normal Form, i.e. with several positive (resp. negative) influences in which the sources are interpreted by a conjunction [FMRS18tcbb].
influence ::= 
inputs ::= 
enumeration ::= 
 _ 
4.1. Influence editor
add_influence(influence).
adds influence rules to the current set of influences.
This command is implicit: influence rules can be added directly in
influence models.
Example.
This example shows a positive influence on b with both a positive source, a, and a negative source (i.e. an inhibitor), b.
The positive and negative sources have opposite effects on the force of the influence,
but do not change the nature of the influence, i.e. the activation of b.
This motivates the use of the positive Boolean semantics which simply ignores the inhibitors and the forces.
On the other hand, the negative Boolean semantics interprets the inhibitors as negative conditions for the influence.
The second influence in this example is added to show the difference between both Boolean dynamics.
biocham:
parameter(k=1).biocham:
present(a,1).biocham:
MA(k)/(1+b) for a / b > b.biocham:
list_ode.[0] d(b)/dt=k*a/(1+b) [1] d(a)/dt=0
biocham:
numerical_simulation.biocham:
plot.biocham:
absent(b). absent(c).biocham:
MA(k)/(1+a) for b / a > c.biocham:
option(boolean_semantics:positive).biocham:
generate_ctl_not.reachable(stable(a)) reachable(stable(b)) reachable(stable(c)) reachable(steady(not c)) reachable(not b) checkpoint2(a,b) checkpoint2(a,c) checkpoint2(b,c) checkpoint2(not c,b) checkpoint2(not b,c)
biocham:
option(boolean_semantics:negative).biocham:
generate_ctl_not.reachable(stable(a)) reachable(stable(b)) reachable(stable(not c)) reachable(not b) checkpoint2(a,b) checkpoint2(not c,b)
list_influences.
lists the current set of influence rules.
influence_model.
creates a new influence model by inferring
the influences between all molecular objects of the current hybrid model
reaction_model.
creates a new reaction model by inferring
the reactions between all molecular objects of the current hybrid model
4.2. Influence graph
Biocham can manipulate three notions of influence graph of a model.
First, the influence hypergraph of an influence network or a reaction network depicts the structure of the influences between the constituants of a model. It is represented by a bipartite moleculeinfluence graph.
Second, the influence graph of a reaction or influence network is a signed directed simple graph between the molecular species.
This graph an abstraction defined from the stoichiometry of the reactions, which is equivalent, under general conditions, to the influence graph defined by the signs of the Jacobian matrix of the ODEs [FS08fmsb,FGS15tcs].
Third, the multistability graph is a multigraph variant of the influence graph in which the influence arcs are labelled by the reactions from which they originate.
This labelled graph can be used for checking very efficiently necessary conditions for the existence of oscillations and multiple steady states [BFS18jtb] (
check_multistability/0
).influence_hypergraph.
builds the hypergraph of the influence rules of the current model.
draw_influence_hypergraph.
Draws the hypergraph of the influence rules of the current model. The hypergraph distinguishes each influence. Equivalent to
influence_hypergraph. draw_graph
.
influence_graph.
builds the influence graph between molecular species without distinguishing the different influences.
draw_influences.
Draws the influence graph between molecular species of the current model, by merging the different influences. Equivalent to
influence_graph. draw_graph
.
Example.
biocham:
load(library:examples/mapk/mapk).biocham:
draw_influences.export_lemon_graph(output_file).
exports the current influence or multistability graph to
output_file
in Lemon
graph format (http://lemon.cs.elte.hu/trac/lemon) (adding .lgf
extension if needed). Computes the conservation laws of the model (by search_conservations/0
) in order
to do so.
multistability_graph.
Creates the influence graph of the current model as by
influence_graph/1
but with the arcs labelled with the reactions they originate from. This is used for command check_multistability/0
.
Options.
 force_graph: yesno
 Force the creation of the graph
option(force_graph:no).
Chapter 5
Events
Events can be used to change some parameter values once a condition gets satisfied. This is useful to implement discrete events in a variety of situations.
Events can be used in modelling, for instance for dividing the cell mass by 2 at each division in a model of the cell cycle,
or for creating hybrid automata model, which chain different continuous semantics (ODEs).
Events can also be intensively used to implement stochastic simulators defined by events, and hybrid differentialstochastic simulators [CFJS15tomacs].
add_event(condition, parameter_name_{1} = arithmetic_expression_{1},
...,
parameter_name_{n} =
arithmetic_expression_{n}).
sets up an event that will be fired each time the condition given as first
argument goes from false to true.
This command is effective in numerical simulations only.
Upon firing, the parameters receive new values
computed from the expression.
The initial values of the parameters are restored after the simulation.
Example.
biocham:
'MA'(k)for a=>b.biocham:
parameter(k=1).biocham:
add_event(b>0.5,k=0).biocham:
present(a).biocham:
numerical_simulation(time:2).biocham:
plot.biocham:
numerical_simulation(time:2,method:ssa).biocham:
plot.
Note that the continuous
numerical_simulation/0
engine will
attempt to interpolate linear event conditions as per [doi.org/10.1007/3540453512_19].
add_time_event(condition, parameter_name_{1} = arithmetic_expression_{1},
...,
parameter_name_{n} =
arithmetic_expression_{n}).
Introduce a special kind of event when Time is the only variable, they may be used
for efficiency reason during numerical integration. The only allowed syntax is of the
form: Time > a or Time >= b.
list_events.
lists all the declared events.
Chapter 6
Importing and Exporting BIOCHAM Models
6.1. Biocham files
load(input_file).
acts as the corresponding
load_biocham/1
⁄ load_sbml/1
⁄ load_ode/1
⁄ load_table/1
,
depending on the file extension
(respectively .bc
, .xml
, .ode
, .csv
– assuming no extension is .bc
).
add(input_file).
acts as the corresponding
add_biocham/1
⁄ add_sbml/1
⁄ add_ode/1
⁄ add_table/1
,
depending on the file extension
(respectively .bc
, .xml
, .ode
, .csv
– assuming no extension is .bc
).
load_biocham(input_file).
opens a new model, loads the reaction rules and executes the commands
(with the file directory as current directory)
contained in the given Biocham
.bc
file.
The suffix .bc
is automatically added to the name if such a
file exists.
add_biocham(input_file).
the rules of the given
.bc
file are loaded and
added to the current set of rules.
The commands contained in the file are executed
(with the file directory as current directory).
export_biocham(output_file).
exports the current model into a
.bc
file.
new_model.
opens a new fresh model.
clear_model.
clears the current model.
list_models.
lists all open models.
list_current_models.
lists current models.
list_model.
lists the current Biocham model: reactions, influences, initial state, parameters, events, functions, LTL patterns, options.
select_model({ref_{1}, ...,
ref_{n}}).
selects some models.
set_model_name(name).
changes the current model name.
delete([range_{1}],
...,
[range_{n}]).
deletes the listed elements from the model.
inherits(ref_{1},
...,
ref_{n}).
makes the current model inherit from the given ancestor models.
parametrize.
Replace numeric constants in a model by parameters k1, k2,... initialized with their value.
6.2. SBML and SBMLqual files
load_sbml(input_file).
acts as
load_biocham/1
but importing reactions, parameters and
initial state (and only that!) from an SBML .xml file.
add_sbml(input_file).
acts as
add_biocham/1
but importing reactions, parameters and
initial state (and only that!) from an SBML .xml file.
option(output_to_library:no).
download_curated_biomodel(Id: integer).
Downloads to the current directory the SBML file for the corresponding
(curated) biomodel with numeric ID
Id
(e.g. "5").
Options.
 output_to_library: yesno
 outputs to the library⁄biomodels subfolder
export_sbml(output_file).
exports the current model into an SBML
.xml
file.
load_qual(input_file).
acts as
load_biocham/1
but importing influences and
initial state (and only that!) from an SBML3qual .sbml file.
add_qual(input_file).
acts as
add_biocham/1
but importing influences and
initial state (and only that!) from an SBML3qual .sbml file.
load_ginml(input_file).
acts as
load_biocham/1
but importing influences
from a GINsim .ginml file.
add_ginml(input_file).
acts as
add_biocham/1
but importing influences
from a GINsim .ginml file.
6.3. Ordinary Differential Equations
Biocham can also manipulate ODE systems, import and export ODE systems in XPP syntax, export in LaTeX, and also infer an equivalent reaction network or influence network from the ODEs [FGS15tcs]. This is useful for importing MatLab models in SBML, and using Biocham analyzers on ODE models.
Remark. XPP format has restrictions on names and does not distinguish between upper case and lower case letters. As a consequence, some complex names may be not properly exported in ode files, and when an ode file is imported, the names are read in lower case.
ode ::= 
oderef ::= 
Listing ODEs
list_ode.
returns the set of ordinary differential equations
and initial concentrations (one line per molecule)
associated to the current model.
Example.
biocham:
a=>b.biocham:
list_ode.[0] d(b)/dt=a [1] d(a)/dt= a
Exporting ODEs
export_ode(output_file).
exports the ODE system of the current model.
export_ode_as_latex(output_file).
exports the ODE system of the current model as LaTeX code.
Inferring reaction or influence networks from an ODE file
Example.
biocham:
parameter(k=10).biocham:
MA(k) for 2*a + b => 3*c.biocham:
list_ode.[0] d(c)/dt=3*k*a^2*b [1] d(b)/dt=  (k*a^2*b) [2] d(a)/dt=  (2*k*a^2*b)
biocham:
export_ode('test2.ode').biocham:
load_reactions_from_ode('test2.ode').biocham:
list_model.k*a^2*b for 2*a+b=>3*c. present(c,0). present(b,0). present(a,0). parameter( k = 10 ).
biocham:
load_influences_from_ode('test2.ode').biocham:
list_model.3*(k*a^2*b) for b,a > c. 2*(k*a^2*b) for b,a < a. 1*(k*a^2*b) for b,a < b. present(c,0). present(b,0). present(a,0). parameter( k = 10 ).
biocham:
list_ode.[0] d(b)/dt=  (k*a^2*b) [1] d(a)/dt=  (2*k*a^2*b) [2] d(c)/dt=3*k*a^2*b
option(import_reactions_with_inhibitors:no).
load_reactions_from_ode(input_file).
infers a set of reactions equivalent to an ODE system, and loads it as
load_biocham/1
.
Options.
 import_reactions_with_inhibitors: yesno
 Add inhibitors when inferring reactions.
add_reactions_from_ode(input_file).
infers a set of reactions equivalent to an ODE system, and adds it to the current model as
add_biocham/1
.
Options.
 import_reactions_with_inhibitors: yesno
 Add inhibitors when inferring reactions.
load_influences_from_ode(input_file).
infers a set of influences equivalent to an ODE system, and loads it as
load_biocham/1
.
add_influences_from_ode(input_file).
infers a set of influences equivalent to an ODE system, and adds it to the current model as
add_biocham/1
.
ODE systems
ode_system.
builds the set of ODEs of of the current model.
import_ode(input_file).
imports a set of ODEs.
new_ode_system.
creates an ODE system.
delete_ode_system(name).
deletes an ODE system.
set_ode_system_name(name).
sets the name of the current ODE system.
list_ode_systems.
lists the ODE systems of the current model.
select_ode_system(name).
selects an ODE system
add_ode(ode_{1},
...,
ode_{n}).
If there is a current ODE system, adds the given set of ODEs to it.
If there is no current ODE system,
infers reactions equivalent to the given set of ODEs.
delete_ode(oderef_{1},
...,
oderef_{n}).
removes the given set of ODEs from the current ODE system.
init(name_{1} = arithmetic_expression_{1},
...,
name_{n} =
arithmetic_expression_{n}).
sets the initial value of a variable in the current set of ODEs.
add_reactions_from_ode_system.
adds a set of reactions equivalent to the current ODE system.
Options.
 import_reactions_with_inhibitors: yesno
 Add inhibitors when inferring reactions.
add_influences_from_ode_system.
adds a set of influences equivalent to the current ODE system.
load_reactions_from_ode_system.
Replaces the current reactions with those from
add_reactions_from_ode_system/0
.
Options.
 import_reactions_with_inhibitors: yesno
 Add inhibitors when inferring reactions.
load_influences_from_ode_system.
Replaces the current reactions with those from
add_influences_from_ode_system/0
.
Part II: Qualitative Analysis and Synthesis
Chapter 7
Static Analyses
7.1. Dimensional analysis
Dimensional analysis infers dimensions for model parameters and checks their consistency.
In Biocham, only time and volume dimensions are considered, without units.
The dimension of a molecular concentration is volume$^{1}$.
The dimension of a kinetic expression (i.e. reaction rate or influence force) is time$^{1}$.
Warning: numbers have no dimension so kinetic parameters should be used in kinetic expressions instead of directly numbers.
list_dimensions.
Infers the time and volume dimensions of all parameters.
Example.
biocham:
MM(v, k) for A => B.biocham:
parameter(k = 2, v = 3).biocham:
list_dimensions.v has dimension time^(1).volume^(0) k has dimension time^(0).volume^(1)
set_dimension(P: parameter_name, U: number, V: number).
Declare dimension time^
U
. volume^V
for parameter P
clear_dimensions.
Clear all declared dimensions
clear_dimension(P: parameter_name).
Clear declared dimension for parameter
P
7.2. Conservation laws and invariants
Petri Net Placeinvariants provide linear conservation laws for the differential semantics of a reaction network.
They can be verified or computed with the commands below. Note that such invariants are useful information but it may be not a good idea to use them to reduce
the dimensionality of the system since it may introduce numerical instabilities.
add_conservation(Conservation: solution).
declares a new mass conservation law for all molecules given with the
corresponding weight in
Conservation
. During a numerical
simulation, one of those variables will be
eliminated thanks to this conservation law. Be careful if you declare
conservation laws and then plot the result of a previous simulation, the
caption might not be correct.
When added, the conservation law will be checked against the reactions
(i.e. purely from stoichiometry), if that fails against the kinetics.
Since these checks are not complete, even a failure will be accepted with
a warning.
delete_conservation(Conservation: solution).
removes the given mass conservation law.
delete_conservations.
removes all mass conservation laws.
list_conservations.
prints out all the mass conservation laws.
check_conservations.
checks all conservation laws against reactions, and if necessary kinetics
(see also
add_conservation/1
).
option(conservation_size:4).
search_conservations.
computes and displays the Pinvariants of the network up to the maximal
size given by the option conservation_size and defaulting to 4. Such
Pinvariants are particular mass conservation laws that are independent
from the kinetics.
Options.
 conservation_size: integer
 Set the maximal stoichiometric size for a Pinvariant.
7.3. Detecting model reductions
This section describes commands to detect model reduction
relationships between reaction models based solely on the
structure of their reaction graph [GSF10bi].
The commands below check the existence of a subgraph
epimorphism (SEPI) [GFS14dam], i.e. a graph reduction
from one graph to a second graph, obtained by deleting and⁄or
merging species and⁄or reactions. A SAT solver is used to
solve this NPcomplete problem.
extremal_sepi ::= 
 Search standard reduction.no  Search reduction that minimises bottom (i.e., the number of deletions).minimal_deletion  Search reduction that maximises bottom (i.e., the number of deletions).maximal_deletion 
merge_restriction ::= 
 Search standard reduction.no  Search reduction with local twoneighbour restriction (cf. report chapter 5).neighbours  Search reduction with old twoneighbour restriction (prior to 2019).old 
option(mapping_restriction:[]).
option(merge_restriction:no).
option(timeout:180).
option(all_reductions:no).
option(distinct_species:no).
option(max_nb_reductions:200).
option(extremal_sepi:no).
option(stats:no).
search_reduction(FileName1: input_file, FileName2: input_file).
checks whether there exists one model reduction from a
first Biocham model to a second model, and returns the
first model reduction found, as a description of a graph
morphism (SEPI) from the reaction graph of the first model
to the second. Optionally, a partial mapping of the form
['label1' > 'label2', 'label3' > 'deleted'] can
be given to restrict the search to reductions satisfying
this mapping of molecular names.
Another option restricts the merge operation to nodes
sharing at least one neigbour in the reaction graph.
A timeout option can be given in seconds (default is 180s).
Note that the files FileName1 and FileName2 will be loaded,
therefore any previous model will be overwritten and some
biocham commands might be executed.
Options.
 mapping_restriction: [basic_influence_{1}, ..., basic_influence_{n}]
 enforce the given mappings in the SEPI
 merge_restriction: merge_restriction
 restrict merges according to some criterion
 timeout: number
 timeout for the (Max)SAT solver
 all_reductions: yesno
 specifies if solver is looking for all SEPI reductions or not
 distinct_species: yesno
 specifies if solver is listing only SEPIs with distinct species (cf. report section 6.6)
 max_nb_reductions: number
 limits the number of SEPI reductions the solver is looking for
 extremal_sepi: extremal_sepi
 defines the type of reduction searched (cf. report chapter 3)
 stats: yesno
 display computation time
Example.
Detecting MichaelisMenten reductions:
biocham:
load(library:examples/sepi/MM1.bc).biocham:
list_model.MA(1) for E+S=>ES. MA(1) for ES=>E+S. MA(1) for ES=>E+P.
biocham:
load(library:examples/sepi/MM2.bc).biocham:
list_model.MA(1) for B+A=>C+A.
biocham:
search_reduction(library:examples/sepi/MM2.bc, library:examples/sepi/MM1.bc).no sepi found
biocham:
search_reduction(library:examples/sepi/MM1.bc, library:examples/sepi/MM2.bc, mapping_restriction: [E>A]).sepi E > A S > B ES > deleted P > C reaction0 > reaction0 reaction1 > deleted reaction2 > reaction0 Number of reductions: 1
7.4. Pattern reduction
This section describes commands to rewrite a reaction graph by reducing some specified patterns (for example MichaelisMenten patterns). The idea behind this function is to use graph rewriting before searching SEPIs (with the command
search_reduction/2
. It is also possible to expand reduced MichaelisMenten patterns, in order to exhibit the intermediary step of reversible complexation enzymesubstrate.option(michaelis_menten:yes).
option(r_1:yes).
option(r_2:no).
option(ep:no).
option(enzyme:yes).
option(hill_reaction:no).
option(partial_hill_reaction:no).
option(double_michaelis_menten:no).
option(michaelis_menten_expansion:no).
yesnomaybe ::= 
 Value yes.yes  Value no.no  Value maybe.maybe 
pattern_reduction(Input_file: input_file, Output_file: input_file).
The patterns to be searched into the input graph (the graph associated to the Biocham model given as
Input_file
) are specified by the options. For a pattern, being found into a graph means that there is a subgraph isomorphism from the graph to the pattern and that only the vertices sent to the image vertices enzyme, substrate, product can have edges with other vertices in the graph (for example it constrains the complex not to intervene in any other reaction). The number of found patterns is printed in the terminal. In the output model (written in a new file .bc, with name specified by the user as Output_file
), all occurrences of these patterns are reduced to the onereaction MichaelisMenten pattern, with a doublearrow from the enzyme. The michaelis_menten_expansion option has the opposite effect to expand reduced MichaelisMenten patterns
Options.
 michaelis_menten: yesno
 specifies if reducing MichaelisMenten patterns
 r_1: yesnomaybe
 R1 reaction (from the commplex ES to enzyme E and substrate S) in the searched MichaelisMenten pattern
 r_2: yesnomaybe
 R2 reaction (from the product P and enzyme E to the complex EP or ES) in the searched MichaelisMenten pattern
 ep: yesnomaybe
 irreversible reaction from the complex ES to the complex EP in the searched MichaelisMenten pattern
 hill_reaction: yesno
 specifies if reducing Hill patterns
 partial_hill_reaction: yesno
 specifies if reducing partial Hill patterns
 double_michaelis_menten: yesno
 specifies if reducing double MichaelisMenten
 enzyme: yesno
 specifies if the reduced MichaelisMenten pattern (that will appear in the output graph) is with enzyme or no
 michaelis_menten_expansion: yesno
 expansion of MichaelisMenten patterns (form the reduced form with enzyme to form with complex ES and reaction R1), should be used with option michaelis_menten:no to avoid confusion
Example.
Example of reduction of a MichaelisMenten pattern with reactions R1, R2 and complex EP.
Input graph :
biocham:
load(library:examples/sepi/MM_variants/MM_ESP.bc).biocham:
option(michaelis_menten:yes).biocham:
option(r_1:maybe,r_2:maybe,ep:maybe).biocham:
pattern_reduction(library:examples/sepi/MM_variants/MM_ESP.bc, MM_ESP_reduced.bc).Number of Michaelis Menten patterns: 1
biocham:
load(MM_ESP_reduced.bc).
Output graph (reduced MichaelisMenten pattern) :
Example.
Hill pattern (can be reduced with option hill_reaction) :
Example.
Partial Hill pattern (can be reduced with option partial_hill_reaction) :
Example.
Double MichaelisMenten pattern (can be reduced with option double_michaelis_menten) :
Example.
Example of expansion of reduced MichaelisMenten in a MAPK cascade.
biocham:
load(library:examples/sepi/mapk0.bc).biocham:
draw_reactions.biocham:
option(michaelis_menten:no).biocham:
option(michaelis_menten_expansion:yes).biocham:
pattern_reduction(library:examples/sepi/mapk0.bc, mapk0_expanded.bc).Number of Michaelis Menten patterns: 2
biocham:
load(mapk0_expanded.bc).biocham:
draw_reactions.7.5. Tropical algebra equilibrations
Tropical algebra can be used to reason about orders of magnitude of molecular concentrations, kinetic parameters and reactions rates.
While steady state analysis consists in finding the roots of the differential functions associated to all or some molecular species in a CRN,
tropical equilibration consists in finding when the positive and negative preponderant terms of the differential functions have the same order of magnitude
(i.e. same integer logarithm).
The solutions to tropical equilibration problems provide candidates for regimes exhibiting fastslow dynamics leading to model reductions
based on quasisteady states or reactions in quasiequilibrium [SFR14amb].
option(tropical_epsilon:0.1).
option(tropical_max_degree:3).
option(tropical_ignore:{}).
option(tropical_denominator:0).
option(tropical_single_solution:no).
tropicalize.
Try to solve a tropical equilibration problem.
Example.
biocham:
load(library:examples/cell_cycle/Tyson_1991.bc).biocham:
tropicalize(tropical_max_degree: 2).Cdc2+Cdc2Cyclin~{p1,p2}+Cdc2Cyclin~{p1}+Cdc2~{p1} 1 complex invariant(s) found a complete equilibration leading to the rescaling: Cdc2' = ε^( 2) * Cdc2 Cdc2Cyclin~{p1,p2}' = Cdc2Cyclin~{p1,p2} Cdc2Cyclin~{p1}' = ε^( 2) * Cdc2Cyclin~{p1} Cdc2~{p1}' = Cdc2~{p1} Cyclin' = ε^( 4) * Cyclin Cyclin~{p1}' = ε^( 2) * Cyclin~{p1} found a complete equilibration leading to the rescaling: Cdc2' = ε^( 3) * Cdc2 Cdc2Cyclin~{p1,p2}' = Cdc2Cyclin~{p1,p2} Cdc2Cyclin~{p1}' = ε^( 2) * Cdc2Cyclin~{p1} Cdc2~{p1}' = ε^( 1) * Cdc2~{p1} Cyclin' = ε^( 3) * Cyclin Cyclin~{p1}' = ε^( 2) * Cyclin~{p1} found a complete equilibration leading to the rescaling: Cdc2' = ε^( 4) * Cdc2 Cdc2Cyclin~{p1,p2}' = Cdc2Cyclin~{p1,p2} Cdc2Cyclin~{p1}' = ε^( 2) * Cdc2Cyclin~{p1} Cdc2~{p1}' = ε^( 2) * Cdc2~{p1} Cyclin' = ε^( 2) * Cyclin Cyclin~{p1}' = ε^( 2) * Cyclin~{p1} No (more) complete equilibration
Options.
 tropical_epsilon: number
 Used for computing degrees (must be < 1)
 tropical_max_degree: number
 Max absolute value of the degree in epsilon as a power of 2
 tropical_ignore: {object_{1}, ..., object_{n}}
 set of species to ignore for equilibration ({} to ignore none, {1} to ignore all unbalanced)
 tropical_denominator: number
 Look for solutions that are integers⁄(tropical_denominator+1)
 tropical_single_solution: yesno
 Stop after finding one solution
 stats: yesno
 display computation time
7.6. Multistability analysis
The existence of oscillations and multiple nondegenerate steady states in the differential dynamics of a reaction or influence model
can be checked very efficiently by checking necessary conditions for those properties
in the multistability graph of the model (see
multistability_graph/0
), namely the existence of respectively negative and positive circuits in that graph [BFS18jtb]. Some more computationally expensive conditions are optional.
option(test_permutations:no).
option(test_transpositions:no).
check_multistability.
Checks the possibility of multistability in the continuous semantics of the current model. This commands checks a necessary condition for the existence of multiple nondegenerate (i.e. with no variable equal to 0) steady states.
Options.
Example.
biocham:
load("library:biomodels/BIOMD0000000040.xml").biocham:
multistability_graph.biocham:
draw_graph(left_to_right:yes).biocham:
check_multistability.There may be nondegenerate multistationarity, positive circuit detected.
biocham:
check_multistability(test_transpositions:yes).There is no multiple steady states with nonzero values, no positive circuit (using permutation of species: [[Br,HBrO2],[Ce,Ce],[HBrO2,Br]]) (using change of sign for species : HBrO2 )
check_oscillations.
Checks a necessary condition for the periodic behaviour of the current model.
Chapter 8
Boolean Dynamics, Verification and Synthesis
8.1. Boolean attractors
The following commands refer to the Boolean dynamics of BIOCHAM models which can be either positive (i.e. without negation, the inhibitors are ignored)
or negative (the inhibitors of a reaction or an influence are interpreted as negative conditions) [FMSR16cmsb].
list_stable_states.
lists stable steady states of the state transition graph corresponding
to the current boolean semantics of the current influence model.
Options.
 boolean_semantics: boolean_semantics
 Use positive or negative boolean semantics for inhibitors.
list_tscc_candidates.
lists possible representative states of Terminal Strongly Connected
Components (TSCC) of the state transition graph corresponding to the
positive semantics of the current influence model.
8.2. Computation Tree Logic (CTL) formulae
The Computation Tree Logic CTL can be used to express qualitative properties of the behavior of a network in a given (set of) initial states, at the boolean level of abstraction [CCDFS04tcs].
CTL extends propositional logic with modal operators to specify where (E: on some path, A: on all paths)
and when (X: next state, F: finally at some future state, G: globally on all states, U: until a second formula is true, R: release)
a formula is true. As in any logic, these modalities can be combined with logical connectives and imbricated,
with the only restriction that a temporal operator must be immediately preceded by a path quantifier,
e.g.
EF(AG(p))
which expresses the reachability of a stable state where protein p
will always remain present.
CTL is well suited to analyze attractors, their reachability and transient properties [TFFT16bi].
It is worth noting however that CTL reachability properties are purely factual and not necessarily causal.The syntax of CTL formulas and some useful abbreviations are defined by the following grammar:
ctl ::= 
reachable(f) is a shorthand for EF(f)
mustreach(f) is a shorthand for AF(f)
steady(f) is a shorthand for EG(f)
stable(f) is a shorthand for AG(f)
checkpoint(f, g) is a shorthand for not EU(not f, g)). Note that such a factual property does not imply any causality relationship from f to g.
checkpoint2(f, g) is a shorthand for EF(g)⁄\checkpoint(f,g)
oscil2(f) is a shorthand for EU(not(f), f ⁄\ EU(f, not(f) ⁄\ EU(not(f), f))) which is true if f has at least two peaks on one path
oscil3(f) is a shorthand for EU(not(f), f ⁄\ EU(f, not(f) ⁄\ EU(not(f), f ⁄\ EU(f, not(f) ⁄\ EU(not(f), f))))) which is true if f has at least three peaks on one path
oscil(f) is a shorthand for oscil3(f) ⁄\ EG(EF(f) ⁄\ EF(not f)) which is a necessary (not sufficient) condition for infinite oscillations in CTL
The qualitative behavior of a network can be specified by a set of CTL formulas using the following commands:
list_ctl.
Prints out all formulae from the current CTL specification.
expand_ctl(Formula: ctl).
shows the expansion in CTL of a formula with patterns.
add_ctl(Formula: ctl).
Adds a CTL formula to the currently declared CTL specification.
delete_ctl(Formula: ctl).
Removes a CTL formula to the currently declared CTL specification.
delete_ctl.
Removes all formulae from the current CTL specification.
The set of CTL formulas of some simple form that are true in the current set of initial states
can also be automatically generated as a specification of the network using the following commands.
generate_ctl(Formula: ctl).
adds to the CTL specification all the CTL formulas that are true, not subsumed by another CTL formula in the specification, and that match the argument formula in which the names that are not molecules are treated as wildcards and replaced by distinct molecules of the network. This command is a variant with subsumption test of
add_ctl
if all names match network molecule names, otherwise it enumerates all m^v instances (where m is the number of molecules and v the number of wildcards in the formula).
generate_ctl.
adds to the CTL specification all reachable stable, reachable steady, reachable, checkpoint2, and oscil formulas on all molecules (using this order of enumeration for performing subsumption checks).
generate_ctl_not.
adds to the CTL specification all reachable stable , reachable stable not, reachable steady, reachable steady not, reachable, reachable not, checkpoint2, and oscil formulas on all molecules.
cleanup_ctl.
cleans up the CTL specification by removing redundant formulae such as reachable(steady(A)) if reachable(stable(A)) is true, etc.
Example.
biocham:
a=>b.biocham:
b=>c.biocham:
c=>d.biocham:
present(b).biocham:
make_absent_not_present.biocham:
generate_ctl(checkpoint2(x,d)).checkpoint2(b,d) checkpoint2(c,d)
biocham:
generate_ctl_not.reachable(stable(d)) reachable(stable(not a)) reachable(stable(not b)) reachable(stable(not c)) reachable(steady(b)) reachable(steady(c)) reachable(steady(not d)) checkpoint2(b,c) checkpoint2(b,d) checkpoint2(c,d) checkpoint2(not a,c) checkpoint2(not d,c) checkpoint2(not a,d) checkpoint2(not c,d) checkpoint2(not a,not b) checkpoint2(not c,not b) checkpoint2(not d,not b) oscil(c)
CTL properties can be efficiently verified with modelcheckers. Biocham uses the NuSMV model checker [nusmv] for which the following options can be specified:
boolean_semantics ::= 
 positive  negative 
The default value is
negative
option(boolean_semantics:negative).
nusmv_initial_states ::= 
 all  some  all_then_some 
The default value is
all
. The value
all_then_some
tries for all initial states, and if that fails, for
someoption(nusmv_initial_states:all).
option(nusmv_counter_example:no).
option(query:current_spec).
check_ctl.
evaluates the current CTL specification (i.e., the conjunction of all
formulae of the current specification) or the content of option
query
on the current model by calling the NuSMV modelchecker.
As is usual in ModelChecking, the query is evaluated for all possible
initial states (Ai
in Biocham v3). This can be changed via
the nusmv_initial_states
option.
Options.
 query: ctl
 Query to evaluate instead of the current CTL specification.
 nusmv_initial_states: nusmv_initial_states
 Consider that a query is true if verified for all⁄some initial states.
 nusmv_counter_example: yesno
 Compute a counterexample for a query when possible (default value is "no").
 boolean_semantics: boolean_semantics
 Use positive or negative boolean semantics for inhibitors.
Example.
biocham:
present(a).biocham:
absent(b).biocham:
a=>b.biocham:
a+b=>a.biocham:
check_ctl(query: EX(not(a) \/ EG(not(b))), nusmv_counter_example:yes).EX(not a\/EG(not b)) is true
biocham:
check_ctl(query: EG(not(b)), nusmv_counter_example:yes).Trace: a b TRUE FALSE EG(not b) is false
biocham:
check_ctl(query:reachable(b),nusmv_counter_example:yes).reachable(b) is true
biocham:
generate_ctl.reachable(stable(b)) reachable(steady(a)) checkpoint2(a,b) oscil(b)
biocham:
check_ctl.current_spec is true
export_nusmv(output_file).
exports the current Biocham set of reactions and initial state in an SMV
.smv
file.
Options.
 boolean_semantics: boolean_semantics
 Use positive or negative boolean semantics for inhibitors.
8.3. Model reduction from CTL specification
This section describes commands to reduce a reaction model by deleting species and reactions, while preserving a CTL specification of the behaviour.
reduce_model(Query: ctl).
Deletes reactions as long as the specification of the behavior given by a CTL formula passed as
argument remains satisfied in the current initial state.
reduce_model.
Same as above using the current CTL
specification of the behavior.
Example.
biocham:
a=>b.biocham:
b=>c.biocham:
c=>d.biocham:
present(b).biocham:
make_absent_not_present.biocham:
generate_ctl.reachable(stable(d)) reachable(steady(b)) reachable(steady(c)) checkpoint2(b,c) checkpoint2(b,d) checkpoint2(c,d) oscil(c)
biocham:
reduce_model.removed [MA(1) for a=>b]
Example.
Reduction of the MAPK model with respect to the output reachability property only:
all dephosphorylation reactions can be removed, resulting in a different bifurcation diagram with full memory effect.
biocham:
load('library:examples/mapk/mapk.bc').biocham:
reduce_model(reachable('MAPK~{p1,p2}')).removed [MA(1) for 'RAFRAFK'=>RAF+RAFK,MA(1) for 'RAF~{p1}'+RAFPH=>'RAF~{p1}RAFPH',MA(1) for 'RAF~{p1}RAFPH'=>'RAF~{p1}'+RAFPH,MA(1) for 'MEKRAF~{p1}'=>MEK+'RAF~{p1}',MA(1) for 'MEK~{p1}RAF~{p1}'=>'MEK~{p1}'+'RAF~{p1}',MA(1) for MEKPH+'MEK~{p1}'=>'MEK~{p1}MEKPH',MA(1) for 'MEK~{p1}MEKPH'=>MEKPH+'MEK~{p1}',MA(1) for MEKPH+'MEK~{p1,p2}'=>'MEK~{p1,p2}MEKPH',MA(1) for 'MEK~{p1,p2}MEKPH'=>MEKPH+'MEK~{p1,p2}',MA(1) for 'MAPKMEK~{p1,p2}'=>MAPK+'MEK~{p1,p2}',MA(1) for 'MAPK~{p1}MEK~{p1,p2}'=>'MAPK~{p1}'+'MEK~{p1,p2}',MA(1) for MAPKPH+'MAPK~{p1}'=>'MAPK~{p1}MAPKPH',MA(1) for 'MAPK~{p1}MAPKPH'=>MAPKPH+'MAPK~{p1}',MA(1) for MAPKPH+'MAPK~{p1,p2}'=>'MAPK~{p1,p2}MAPKPH',MA(1) for 'MAPK~{p1,p2}MAPKPH'=>MAPKPH+'MAPK~{p1,p2}',MA(1) for 'RAF~{p1}RAFPH'=>RAF+RAFPH,MA(1) for 'MEK~{p1}MEKPH'=>MEK+MEKPH,MA(1) for 'MEK~{p1,p2}MEKPH'=>'MEK~{p1}'+MEKPH,MA(1) for 'MAPK~{p1}MAPKPH'=>MAPK+MAPKPH,MA(1) for 'MAPK~{p1,p2}MAPKPH'=>'MAPK~{p1}'+MAPKPH]
biocham:
dose_response('RAFK',1e5,1e3, time:200, show:'MAPK~{p1,p2}').8.4. Model revision from CTL specification
This section describes commands to revise a reaction model in order to satisfy a given CTL specification of the behavior. The revision algorithm is described in [CCFS05tcsb]
revise_model(Query: ctl).
Use theoryrevision on the current model to satisfy the query given as
argument.
revise_model.
Use theoryrevision as above, using the currently defined CTL
specification.
8.5. PAC Learning influence models from Boolean traces
Implementation of Leslie Valiant's PACLearning [doi.org/10.1215/0961754X2872666] algoritm for kCNF formulae [doi.org/10.1145/1968.1972], for learning influence networks from Boolean traces [CFS17cmsb].
Boolean traces can be generated from Boolean simulations of a reaction or influence model, or from the stochastic simulations.
Threshold values can be specified to transform a stochastic trace in a Boolean trace with the following syntax:
transform ::= 
 none 
generate_traces(NInitialStates: integer, time_{Horizon}, FilePrefix: output_file).
Randomly choses
NInitialStates
initial states and for each,
runs a numerical simulation with method: spn
, method:
ssa
or method: sbn
and
time_{Horizon}
as time
option. All results will be saved
in .csv
files with the given FilePrefix
as prefix
and indexed by the run number (unless the prefix is empty).
Options.
pac_learning({input_file_{1}, ...,
input_file_{n}}).
Uses every timeseries
.csv
file in [input_file_{1}, ...,
input_file_{n}]
as source of samples to learn an influence network.
Options.
 cnf_clause_size: integer
 Maximum size of CNF clauses learnt.
option(data_transform:none).
option(cnf_clause_size:3).
option(boolean_simulation:no).
pac_learning(ModelFile: input_file, NInitialStates: integer, TimeHorizon: integer).
Loads
ModelFile
and runs generate_traces/3
with
the provided arguments and then loads the generated files for
pac_learning/1
.
Example.
Learning LotkaVolterra model from Boolean traces and comparison to the hidden original model.
biocham:
pac_learning(library:examples/lotka_volterra/LVi.bc, 200, 10, boolean_simulation: yes).% Maxmimum K used: 1 % minimum number of samples for h=1: 18 % 109 samples (max h ~ 6.055555555555555) P < P % 0 samples (max h ~ 0) % _ > P % 40 samples (max h ~ 2.2222222222222223) P,R < R % 0 samples (max h ~ 0) % _ > R
biocham:
list_model.k1*R*P for R,P < R. k1*R*P for R,P > P. k2*R for R > R. k3*P for P < P. present(R). present(P). parameter( k1 = 2, k2 = 2, k3 = 1 ).
biocham:
pac_learning(library:examples/circadian_cycle/bernot_comet.bc, 100, 10).% Maxmimum K used: 2 % minimum number of samples for h=1: 54 % 232 samples (max h ~ 4.296296296296297) G,PC < G % 27 samples (max h ~ 0.5) / PC,G > G % 421 samples (max h ~ 7.796296296296297) L < L % 137 samples (max h ~ 2.537037037037037) / L > L % 170 samples (max h ~ 3.1481481481481484) L,PC < PC PC / G < PC % 13 samples (max h ~ 0.24074074074074073) G / PC,L > PC
biocham:
list_model._ / L > L. L < L. _ / G,PC > G. G,PC < G. G / PC,L > PC. PC / G < PC. PC,L < PC.
Options.
Part III: Quantitative Analysis and Synthesis
Chapter 9
Numerical Simulations
9.1. ODE and stochastic simulations
Biocham v4 uses the GNU Scientific Library (GSL http://www.gnu.org/software/gsl/) to perform numerical simulations.
The page
http://www.gnu.org/software/gsl/manual/html_node/SteppingFunctions.html#SteppingFunctions
gives a detailed description of the numerical integration methods and options listed below.
The implicit method
bsimp
for stiff systems is the default one.
method ::= 
 Explicit embedded RungeKutta (2, 3) methodrk2  Explicit 4th order (classical) RungeKutta. Error estimation is carried out by the step doubling method. For more efficient estimate of the error, use the embedded methods described belowrk4  Explicit embedded RungeKuttaFehlberg (4, 5) method. This method is a good generalpurpose integratorrkf45  Explicit embedded RungeKutta CashKarp (4, 5) method.rkck  Explicit embedded RungeKutta PrinceDormand (8, 9) method.rk8pd  Implicit Gaussian first order RungeKutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method.rk1imp  Implicit Gaussian second order RungeKutta. Also known as implicit midpoint rule. Error estimation is carried out by the step doubling method.rk2imp  Implicit Gaussian 4th order RungeKutta. Error estimation is carried out by the step doubling method.rk4imp  Default method. Implicit BulirschStoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.bsimp  A variablecoefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit AdamsBashforth (predictor) and implicit AdamsMoulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12.msadams  Perhaps the most robust method but may be slow, a variablecoefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of nonlinear equations. Method order varies dynamically between 1 and 5.msbdf  Stochastic simulation of a ContinuousTime Markov Chain,
defined as per Gillespie's algorithm [Gillespie76jcp]. Note that
there is no conversion of the initial state, besides rounding to the nearest
integer.ssa  Stochastic Petri net simulation. Similar to the SSA
algorithm above but with a discrete⁄logical time.spn  Random Petri net simulation run. All transitions are
equiprobable.pn  Stochastic Boolean net simulation. Similar to the SSA
algorithm above but with a discrete⁄logical time and Boolean states.sbn  Random Boolean net simulation run. All transitions are
equiprobable.bn  Rosenbrock method.rsbk 
time ::= 
filter ::= 
 Does no filtering at all.no_filter  Only keeps the points that are an extremum for some variable.only_extrema 
option(method:bsimp).
option(error_epsilon_absolute:1.0e6).
option(error_epsilon_relative:1.0e6).
option(initial_step_size:1.0e6).
option(maximum_step_size:0.01).
option(minimum_step_size:1.0e5).
option(precision:6).
option(time:20).
option(stochastic_conversion:100).
option(stochastic_thresholding:1000).
option(filter:no_filter).
option(stats:no).
numerical_simulation.
performs a numerical simulation from time 0 up to a given time.
Options.
 time: time
 time horizon of the numerical integration
 stats: yesno
 display computation time
 method: method
 method for the numerical solver
 error_epsilon_absolute: number
 absolute error for the numerical solver
 error_epsilon_relative: number
 relative error for the numerical solver
 initial_step_size: number
 initial step size for the numerical solver
 maximum_step_size: number
 maximum step size for the numerical solver, as a fraction of the total time
 minimum_step_size: number
 minimum step size, as a fraction of the total time (used to trigger events)
 precision: number
 precision for the numerical solver
 stochastic_conversion: number
 Conversion factor used to scale 1 mole to the given number of molecules.
 stochastic_thresholding: number
 Do not write (but still compute) stochastic steps below one fraction of the total time
 filter: filter
 filtering function for the trace
continue(time).
Sets the initial state to the last state of the last simulation trace
and continues for
time
time units.
9.2. Hybrid simulations
This section describes a hybrid numerical simulation
method. Some reactions which occurs more often and which have
a lot of particles reactants will be seen as continuous,
and other reactions which occurs lesser and with fewer reactants
particles will be seen as stochastic. Species which appears in at
least one continous reaction and one stochastic reaction will be
called hybrid species, the total particle count of one species is
calculated and will be plot at the end of the simulation.
The default simulation time is 200, and the convertion rate between particle
and concentration is 1. The simulation method currently used is the
Rosenbrock method.
hybrid_static_simulation(ODEFilename, StochFilename).
This is a hybrid simulation with a static partitioning of the reactions.
The first input file contains the informations concerning the continous
reactions. The second input file contains the informations concerning the
stochastic reactions. Please list out all the species initial value
(by present(_,_)) in input files.
hybrid_static_simulation(ODEFilename, StochFilename, OutFilename, Volume, Time).
hybrid_dynamic_simulation(InputFile, PropTresh, PartTresh).
This is a hybrid simulation with a dynamic partitioning of the reactions.
The input file, a normal .bc file, contains all the infomation needed. The first
coefficient will be used in the propencity treshold. The second coefficient
is for the particle count treshold. If a reaction meets the two conditions
based on these tresholds anytime during the simulation, the reaction is seen
as continuous. Otherwise, it is viewed as stochastic. Therefore, if theses
numbers are high, the simulation is more precise but slower. On the contrary,
the simulation is quicker but less precise if the numbers are low. Please list
out all the species initial value (by present(_,_)) in input files.
hybrid_dynamic_simulation(InputFile, OutFilename, Volume, Time, PropTresh, PartTresh).
9.3. Plotting and exporting the result of simulations
axes ::= 

 x  y  xy 
floatorauto ::= 
 auto 
option(show:{}).
option(logscale:'').
option(against:'Time').
option(xmin:auto).
option(ymin:auto).
option(xmax:auto).
option(ymax:auto).
plot.
plots the current trace. After a simulation, the trace is composed of molecular concentrations and userdefined functions over time.
Example.
biocham:
load(library:examples/mapk/mapk).biocham:
numerical_simulation(method:msbdf).biocham:
plot.biocham:
plot(show: {MAPK~{p1,p2}, MEK~{p1,p2}}).biocham:
numerical_simulation(filter:only_extrema).biocham:
plot.Options.
 logscale: axes
 Apply logscaling to the specified axes.
 show: {object_{1}, ..., object_{n}}
 Restricts the plot to the given set of objects and functions; everything will be plotted if the set is empty
 against: object
 Selects the X axis for the plot, defaulting to Time.
 xmin: floatorauto
 Select the axes for the current plot (auto is overwrite to *)
 ymin: floatorauto
 Select the axes for the current plot (auto is overwrite to *)
 xmax: floatorauto
 Select the axes for the current plot (auto is overwrite to *)
 ymax: floatorauto
 Select the axes for the current plot (auto is overwrite to *)
Example.
biocham:
a>b.biocham:
a<a.biocham:
present(a,10).biocham:
numerical_simulation(time:5,method:ssa).biocham:
plot.biocham:
numerical_simulation(time:5,method:msbdf).biocham:
plot.biocham:
plot(logscale:y).biocham:
plot(show:a,against:b).biocham:
plot(show:a,against:b,logscale:xy).export_plot(FileTemplate: output_file).
saves the current trace into two files:
FileTemplate
.csv
and .plot
.
export_plot_to_png(output_file).
plots the current trace in a PNG file
export_plot_to_canvas(output_file, BaseOutputFile: output_file).
plots the current trace as a canvas element in an HTML file
9.4. Variations and bifurcations
This section describes commands for continuously varying the values of molecular species or parameters.
This can be used for drawing doseresponse diagrams
under the hypothesis that the input dose is a catalyst not affected by the system.
It can also be used for drawing some simple forms of bifurcation diagrams
obtained by just varying forth and back an input.
variation(object, Value1: arithmetic_expression, Value2: arithmetic_expression).
Makes a molecule continuously vary from Value1 to Value2 in the time option duration, under the hypothesis that the object is not affected by the other reactions on this time scale. This is achieved by adding a synthesis reaction with appropriate kinetics, removing any previous variation reaction on the object, and setting its initial concentration to Value1.
Options.
 time: time
 duration of the variation.
logarithmic_variation(object, Value1: arithmetic_expression, Value2: arithmetic_expression).
Makes a logarithmic variation. Value1 must be non zero.
Options.
 time: time
 duration of the variation.
clear_variation(object).
Deletes any variation of the object and restores its initial state.
clear_variations.
Deletes all variations.
Example.
biocham:
option(time:6).biocham:
variation(a,1.0e7,0.1).biocham:
logarithmic_variation(b,1.0e7,0.1).biocham:
numerical_simulation.biocham:
plot.biocham:
plot(logscale:y).biocham:
list_model.variation_a0.016666649999999998*1 for _=>a. variation_bb/6*log(999999.9999999) for b=>2*b. present(a,1.0e7). present(b,1.0e7).
dose_response(Dose: object, Value1: arithmetic_expression, Value2: arithmetic_expression).
Draws a doseresponse diagram by linear variation of the initial concentration of the input object (the dose) from Value1 to Value2, and plotting the output object (the response) in relation to the dose.
A set of output objects can be given to draw their respective responses.
The time option should be greater or equal to the simulation time horizon necessary to approximate the steady state response.
The variation of the dose is performed over a time horizon of 10 fold the given duration.
Options.
logarithmic_dose_response(Dose: object, Value1: arithmetic_expression, Value2: arithmetic_expression).
Draws a similar doseresponse diagram by logarithmic variation and plots it on a logarithmic scale for both the dose and the responses. Value1 must be non zero.
Example.
This example shows the doseresponse diagrams of the MAPK model of Huang and Ferrell in BiomModels revealing the amplifier and analogdigital converter functions
of this network at the second and third level of the cascade.
biocham:
load_sbml(library:biomodels/BIOMD0000000009.xml).biocham:
option(time:100, show:{P_KKK,PP_KK,PP_K}).biocham:
dose_response(E1,1e6,1e4).change_parameter_to_variable(Object: parameter_name).
Transforms a parameter in a variable initialized to the parameter value (and creates a parameter named variation_Object). This transformation is reversible.
restore_parameter(object).
Restores a varying parameter to a fixed parameter.
bifurcations(Dose: object, Value1: arithmetic_expression, Value2: arithmetic_expression).
Draws a doseresponse diagram by linear variation from Value1 to Value2 and then from Value2 to Value1 (to see memory effects) of the input object dose, and plots the output object responses.
The time duration should be greater or equal to the simulation time horizon necessary to get the response at steady state (in both the increasing forward and decreasing backward variations which may be different).
The increasing and decreasing variations of the dose are performed over a time horizon of 10 fold duration each.
No that a real bifurcation diagram draws the stable and unstable zeros of the differential function of the response. Those zeros are not computed here.
logarithmic_bifurcations(Dose: object, Value1: arithmetic_expression, Value2: arithmetic_expression).
Draws a similar doseresponse diagram by logarithmic variation using a logarithmic scale for both the dose and the response. Value1 and Value2 must be non zero.
Example.
This example of the MAPK network shows a memory effect in one single level of two phosphorylation cycles.
The hysteresis corresponds to the existence of two stable states in the rate equation.
biocham:
load_sbml(library:biomodels/BIOMD0000000026.xml).biocham:
dose_response(MAPKK, 0, 100, time:1e4, show:Mpp).biocham:
bifurcations(MAPKK, 0, 100, time:1e4, show:Mpp).Chapter 10
Verification of Temporal Behaviors and Parameter Synthesis
10.1. Numerical data time series
Numerical data time series, produced by biological experiments or by simulations, can be stored in
.csv
files, and loaded.
They can also be modified by the following commands.load_table(input_file).
loads the given
.csv
file as a table.
export_table(output_file).
exports the current table into a
.csv
file.
list_tables.
lists the current set of tables.
select_table(Table: name).
selects
Table
to be the current table.
rename_table(name).
renames the current table.
delete_table(name_{1},
...,
name_{n}).
deletes some tables.
list_rows.
lists the rows of the current table.
list_columns.
lists the column names of the current table.
Columns can be designated by their index or their name.
delete_column(column_{1},
...,
column_{n}).
deletes the given columns from the current table.
rename_column(column, name).
renames the given column of the current table.
delete_row(number_{1},
...,
number_{n}).
deletes the given rows from the current table.
list_last_state.
lists the last state of the current simulation trace.
10.2. Temporal logic FOLTL(Rlin) formulae
FirtOrder Linear Time Logic with linear constraints over the reals, FOLTL(Rlin), can be used to specify semiqualitative semiquantitative constraints on the dynamical behavior of the system, in a much more flexible manner than by curves to fit [RBFS11tcs]. The syntax of FOLTL(Rlin) formulas is given below with some useful abbreviations [FT14book]:
The
foltl_magnitude
option (default 5) is the multiplicative factor used for the strong comparison operators <<
and >>
over positive numbers, A<<B
means 5*A<B
.option(foltl_magnitude: 5).
foltl ::= 
 false  true 
foltl_predicate ::= 
time_value_list ::= 
parameter_between_interval ::= 
Furthermore, some useful abbreviations have been predefined with the following FOLTL(Rlin) formula patterns:
ltl_pattern forall(X, Formula) = not(exists(X, not(Formula))).
ltl_pattern reachable(Formula) = F(Formula).
ltl_pattern steady(Formula) = G(Formula).
Curves.
ltl_pattern curve(Molecule, Time_Value_List) = curve(Molecule, Time_Value_list).
Sequence of events [doi.org/MRMFJ08bi].
ltl_pattern occurs(Formula) = F(Formula).
ltl_pattern excludes(Formula) = G(not(Formula)).
ltl_pattern invariates(Formula) = G(Formula).
ltl_pattern weak_sequence(Formula1, Formula2) = F(Formula1 /\ F(Formula2)).
ltl_pattern exact_sequence(Formula1, Formula2) = F(Formula1 /\ X(Formula2)).
ltl_pattern sequence(Formula1, Formula2) = G(U(Formula1,Formula2)).
ltl_pattern consequence(Formula1, Formula2) = G(Formula1 => F(Formula2)).
ltl_pattern implication(Formula1, Formula2) = G(Formula1 => Formula2).
Thresholds, global extrema, amplitudes [FT14book].
ltl_pattern reached(Molecule, Concentration) = F(Molecule >= Concentration).
ltl_pattern unreached(Molecule, Concentration) = F(Molecule <= Concentration).
ltl_pattern inf_amplitude(Molecule, Amplitude) = exists(Concentration, F(Molecule <= Concentration /\ F(Molecule >= Concentration+Amplitude))).
ltl_pattern sup_amplitude(Molecule, Amplitude) = exists(Concentration, G(Molecule >= Concentration /\ Molecule <= Concentration+Amplitude)).
ltl_pattern first_value(Molecule, Concentration) = (Molecule=Concentration).
ltl_pattern last_value(Molecule, Concentration) = F(G(Molecule=Concentration)).
ltl_pattern maximum(Molecule, Concentration) = (G(Molecule <= Concentration) /\ F(Molecule>=Concentration)).
ltl_pattern minimum(Molecule, Concentration) = (G(Molecule >= Concentration) /\ F(Molecule<=Concentration)).
ltl_pattern maximum(Molecule, Concentration, T) = (G(Molecule <= Concentration) /\ F(Molecule>=Concentration /\ Time=T)).
ltl_pattern minimum(Molecule, Concentration, T) = (G(Molecule >= Concentration) /\ F(Molecule<=Concentration /\ Time=T)).
ltl_pattern amplitude(Molecule, Amplitude) = exists(Minimum, exists(Maximum, maximum(Molecule,Maximum) /\ minimum(Molecule, Minimum) /\ Amplitude=MaximumMinimum)).
Local extrema
ltl_pattern increase(Molecule, Concentration) = ((Molecule<=Concentration) /\ X(Molecule>=Concentration)).
ltl_pattern decrease(Molecule, Concentration) = ((Molecule>=Concentration) /\ X(Molecule<=Concentration)).
ltl_pattern increase(Molecule) = exists(Concentration, increase(Molecule, Concentration)).
ltl_pattern decrease(Molecule) = exists(Concentration, decrease(Molecule, Concentration)).
ltl_pattern local_maximum(Molecule, Concentration) = F(peak(Molecule,Concentration)).
ltl_pattern local_maximum(Molecule, Concentration, T) = F(peak(Molecule,Concentration,T)).
ltl_pattern peak(Molecule, Concentration) = (Molecule<Concentration /\ X(decrease(Molecule, Concentration))).
ltl_pattern peak(Molecule, Concentration, T) = (Molecule<Concentration /\ X(decrease(Molecule, Concentration) /\ Time=T)).
ltl_pattern first_peak(Molecule, Concentration, T) = U(decrease(Molecule), U(increase(Molecule), peak(Molecule, Concentration, T))).
ltl_pattern first_peak(Molecule, Concentration) = U(decrease(Molecule), U(increase(Molecule), peak(Molecule, Concentration))).
ltl_pattern last_peak(Molecule, Concentration, T) = F(peak(Molecule, Concentration, T) /\ X(U(decrease(Molecule), U(increase(Molecule), exists(c, G(Molecule=c)))))).
ltl_pattern last_peak(Molecule, Concentration) = F(peak(Molecule, Concentration /\ X(U(decrease(Molecule), U(increase(Molecule), exists(c, G(Molecule=c))))))).
ltl_pattern successive_peaks(Molecule, C1, T1, C2, T2) = F(U(decrease(Molecule), U(increase(Molecule), peak(Molecule,C1, T1) /\ X(first_peak(Molecule, C2, T2))))).
ltl_pattern successive_peaks(Molecule, C1, C2) = F(U(decrease(Molecule), U(increase(Molecule), peak(Molecule,C1) /\ X(first_peak(Molecule, C2))))).
ltl_pattern first_successive_peaks(Molecule, C1, T1, C2, T2) = U(decrease(Molecule), U(increase(Molecule), peak(Molecule,C1, T1) /\ X(first_peak(Molecule, C2, T2)))).
ltl_pattern last_successive_peaks(Molecule, C1, T1, C2, T2) = F(peak(Molecule,C1, T1) /\ X(U(decrease(Molecule), U(increase(Molecule), peak(Molecule,C2, T2) /\ X(U(decrease(Molecule), U(increase(Molecule), exists(c, G(Molecule=c))))))))).
ltl_pattern local_minimum(Molecule, Concentration) = F(base(Molecule,Concentration)).
ltl_pattern local_minimum(Molecule, Concentration, T) = F(base(Molecule,Concentration,T)).
ltl_pattern base(Molecule, Concentration) = (Molecule>Concentration /\ X(increase(Molecule, Concentration))).
ltl_pattern base(Molecule, Concentration, T) = (Molecule>Concentration /\ X(increase(Molecule, Concentration) /\ Time=T)).
ltl_pattern first_base(Molecule, Concentration, T) = U(increase(Molecule), U(decrease(Molecule), base(Molecule, Concentration, T))).
ltl_pattern first_base(Molecule, Concentration) = U(increase(Molecule), U(decrease(Molecule), base(Molecule, Concentration))).
ltl_pattern last_base(Molecule, Concentration, T) = F(base(Molecule, Concentration, T) /\ X(U(increase(Molecule), U(decrease(Molecule), exists(c, G(Molecule=c)))))).
ltl_pattern last_base(Molecule, Concentration) = F(base(Molecule, Concentration /\ X(U(increase(Molecule), U(decrease(Molecule), exists(c, G(Molecule=c))))))).
ltl_pattern successive_bases(Molecule, C1, T1, C2, T2) = F(U(increase(Molecule), U(decrease(Molecule), base(Molecule,C1, T1) /\ X(first_base(Molecule, C2, T2))))).
ltl_pattern successive_bases(Molecule, C1, C2) = F(U(increase(Molecule), U(decrease(Molecule), base(Molecule,C1) /\ X(first_base(Molecule, C2))))).
ltl_pattern first_successive_bases(Molecule, C1, T1, C2, T2) = U(increase(Molecule), U(decrease(Molecule), base(Molecule,C1, T1) /\ X(first_base(Molecule, C2, T2)))).
ltl_pattern last_successive_bases(Molecule, C1, T1, C2, T2) = F(base(Molecule,C1, T1) /\ X(U(increase(Molecule), U(decrease(Molecule), base(Molecule,C2, T2) /\ X(U(increase(Molecule), U(decrease(Molecule), exists(c, G(Molecule=c))))))))).
Periods and delays.
ltl_pattern period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, successive_peaks(Molecule, C1, T1, C2, T2) /\ Period=T2T1)))).
ltl_pattern first_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, first_successive_peaks(Molecule, C1, T1, C2, T2) /\ Period=T2T1)))).
ltl_pattern last_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, last_successive_peaks(Molecule, C1, T1, C2, T2) /\ Period=T2T1)))).
ltl_pattern delay(Molecule1, Molecule2, Delay) = exists(T1, exists(T2, exists(C1, exists(C2, F(peak(Molecule1,C1,T1) /\ first_peak(Molecule2,C2,T2) /\ Delay=T2T1))))).
ltl_pattern first_delay(Molecule1, Molecule2, Delay) = exists(T1, exists(T2, exists(C1, exists(C2, U(decrease(Molecule1), U(increase(Molecule1), peak(Molecule1, C1, T1) /\ first_peak(Molecule2,C2,T2) /\ Delay=T2T1)))))).
ltl_pattern last_delay(Molecule1, Molecule2, Delay) = exists(T1, exists(T2, exists(C1, exists(C2, last_peak(Molecule1,C1,T1) /\ last_peak(Molecule2,C2,T2) /\ Delay=T2T1)))).
ltl_pattern base_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, successive_bases(Molecule, C1, T1, C2, T2) /\ Period=T2T1)))).
ltl_pattern first_base_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, first_successive_bases(Molecule, C1, T1, C2, T2) /\ Period=T2T1)))).
ltl_pattern last_base_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, last_successive_bases(Molecule, C1, T1, C2, T2) /\ Period=T2T1)))).
ltl_pattern base_delay(Molecule1, Molecule2, Delay) = exists(T1, exists(T2, exists(C1, exists(C2, F(base(Molecule1,C1,T1) /\ first_base(Molecule2,C2,T2) /\ Delay=T2T1))))).
ltl_pattern first_base_delay(Molecule1, Molecule2, Delay) = exists(T1, exists(T2, exists(C1, exists(C2, U(increase(Molecule1), U(decrease(Molecule1), base(Molecule1, C1, T1) /\ first_base(Molecule2,C2,T2) /\ Delay=T2T1)))))).
ltl_pattern last_base_delay(Molecule1, Molecule2, Delay) = exists(T1, exists(T2, exists(C1, exists(C2, last_base(Molecule1,C1,T1) /\ last_base(Molecule2,C2,T2) /\ Delay=T2T1)))).
FOLTL(Rlin) formulas are evaluated on traces, i.e. numerical data time series either generated by simulation or loaded from biological experiments. FOLTL(Rlin) formulas may contain free variables, in which case they are called constraints. While a closed formula (i.e. without free variable) is either true or false on a trace, a formula with free variable is either satisfiable (i.e. true for some valuations of the variables) or unsatisfiable (i.e. false for any valuation).
validity_domain(Formula: foltl).
solves a FOLTL(Rlin) constraint on the current trace, i.e. computes
the validity domain for the free variables
that make the formula true on the numerical trace.
The formula is false if the validity domain of one of its free variables is empty, and satisfiable otherwise.
Example.
biocham:
a => b.biocham:
present(a,10).biocham:
numerical_simulation. plot.biocham:
validity_domain(G(Time < T => a > b)).T<=1.48828
Options.
 foltl_magnitude: number
 order of magnitude for << and >> operators
satisfaction_degree(Formula: foltl, [variable_objective_{1}, ...,
variable_objective_{n}]).
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).
Options.
 foltl_magnitude: number
 order of magnitude for << and >> operators
ltl_pattern(function_prototype_{1} = foltl_{1},
...,
function_prototype_{n} =
foltl_{n}).
sets the definition of patterns.
list_ltl_patterns.
lists all known LTL patterns.
expand_ltl(Formula: foltl).
shows the expansion in LTL of a formula with patterns.
delete_ltl_pattern(functor_{1},
...,
functor_{n}).
deletes some LTL patterns. Either arity is given, or all LTL patterns with
the given functor are deleted.
10.3. Parameter sensitivity, robustness and parameter optimization w.r.t. FOLTL(Rlin) properties
The continuous satisfaction degree of an FOLTL(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 FOLTL(Rlin) formula [RBFS11tcs] or even maximizing the margins and the robustness [FS18cmsb].
option(robustness_samples: 100).
option(robustness_relative_error: 0.05).
option(robustness_coeff_var: 0.1).
robustness(Formula: foltl, [parameter_name_{1}, ...,
parameter_name_{n}], [variable_objective_{1}, ...,
variable_objective_{n}]).
computes the robustness degree as defined in [RBFS11tcs], with
respect to formula
Formula
, for the list of parameters
[parameter_name_{1}, ...,
parameter_name_{n}]
and with list of objectives for the free variables
of Formula
given in [variable_objective_{1}, ...,
variable_objective_{n}]
.
The robustness is the average satisfaction degree (truncated to one) estimated by sampling of a normal perturbation law with
a coefficient of variation (stddev⁄mean) given by the corresponding option,
and centered around the current parameter values.
Example.
biocham:
k for a => b. present(a,10). parameter(k=0.7).biocham:
robustness(G(Time < T => a > b), [k], [T > 5]).Time: 1.747 s Robustness degree: 1.0
Options.
 time: time
 time horizon of the numerical integration
 method: method
 method for the numerical solver
 robustness_coeff_var: number
 coefficient of variation of the normal law
 robustness_samples: integer
 number of samples used for averaging
 robustness_relative_error: number
 relative sampling error used to stop estimation of the robustness
 error_epsilon_absolute: number
 absolute error for the numerical solver
 error_epsilon_relative: number
 relative error for the numerical solver
 initial_step_size: number
 initial step size for the numerical solver
 maximum_step_size: number
 maximum step size for the numerical solver, as a fraction of time
 precision: number
 precision for the numerical solver
option(cmaes_init_center: no).
option(cmaes_log_normal: no).
option(cmaes_stop_fitness: 0.0001).
search_parameters(Formula: foltl, [parameter_between_interval_{1}, ...,
parameter_between_interval_{n}], [variable_objective_{1}, ...,
variable_objective_{n}]).
tries to satisfy a FOLTL(Rlin) constraint by varying the parameters listed
in
[parameter_between_interval_{1}, ...,
parameter_between_interval_{n}]
.
This search command uses the stochastic optimization procedure CMAES
(covariance matrix adaptation evolution strategy)
with the continuous satisfaction degree (truncated to 1 or not according to stop option) of the given FOLTL(Rlin) property as
fitness function.
Example.
biocham:
k for a => b. present(a,10).biocham:
search_parameters(G(Time < T => a > b), [0 <= k <= 2], [T > 5]).Time: 1.687 s Stopping reason: Fitness: function value 9.74e01 <= stopFitness (1.00e04) Best satisfaction degree: 38.488281 [0] parameter(k=0.11795528701204568)
biocham:
search_parameters(G(Time < T => a > b), [0 <= k <= 2], [T > 5]).Time: 1.63 s Stopping reason: Fitness: function value 5.98e01 <= stopFitness (1.00e04) Best satisfaction degree: 2.488281 [0] parameter(k=0.8243647454738888)
Options.
 time: time
 time horizon of the numerical integration
 method: method
 method for the numerical solver
 error_epsilon_absolute: number
 absolute error for the numerical solver
 error_epsilon_relative: number
 relative error for the numerical solver
 initial_step_size: number
 initial step size for the numerical solver
 maximum_step_size: number
 maximum step size for the numerical solver, as a fraction of time
 precision: number
 precision for the numerical solver
 cmaes_init_center: yesno
 initialize parameter values at the center of their domain instead of their current value
 cmaes_log_normal: yesno
 search on the log of the parameters (better if very different orders of magnitude)
 cmaes_stop_fitness: number
 stop when distance to the objective is less than this value (use 0 by default and negative values in [1,0] to maximize margins)
search_parameters([parameter_between_interval_{1}, ...,
parameter_between_interval_{n}], Conditions).
similar to
search_parameters/3
but uses the list of
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.
Example.
biocham:
ka for _ => a.biocham:
kb for _ => b.biocham:
parameter(ka=1, kb=1).biocham:
search_parameters([0 <= ka <= 3], [(G(ab<x/\ba<x), [], [x > 10]), (G(ab<y/\ba<y), [kb=2], [y > 10])]).Time: 5.459 s Stopping reason: Fitness: function value 4.40e05 <= stopFitness (1.00e04) Best satisfaction degree: 0.999956 [0] parameter(ka=1.5000022018101784)
Options.
 time: time
 time horizon of the numerical integration
 method: method
 method for the numerical solver
 error_epsilon_absolute: number
 absolute error for the numerical solver
 error_epsilon_relative: number
 relative error for the numerical solver
 initial_step_size: number
 initial step size for the numerical solver
 maximum_step_size: number
 maximum step size for the numerical solver, as a fraction of time
 precision: number
 precision for the numerical solver
 cmaes_init_center: yesno
 initialize parameter values at the center of their domain instead of their current value
 cmaes_log_normal: yesno
 search on the log of the parameters (better if very different orders of magnitude)
 cmaes_stop_fitness: number
 stop when distance to the objective is less than this value (use negative to enforce robustness)
Chapter 11
Synthesis of Reaction Networks
11.1. Synthesis from mathematical expressions and simple programs
In this simple program syntax each statement must be terminated by ',' or '.'.
statement ::= 
assignment ::= 
variable ::= 
expression ::=  The support of arithmetic_expression is not complete, but most of them are supported. 
Statements terminated by ',' will execute parallelly with next statement.
lhs ::= 
rhs ::= 
control_flow_keyword ::= 
 if_tag  When using these two keywords, the rhs serves as the condition.while  else_tag  endif  When using these three keywords, the rhs does not matter.endwhile 
The operators 'pre' and 'post' are used to describe the procedural structure of a program,
so if not specified, the list element will execute parallelly, which makes it possible to describe sequential and parallel part of a program.
To mix the sequential and parallel computation, statements which executed parallelly need to add 'para' before the expression, except the last statement.
Currently it is user's responsibility to ensure the execution of last statement finishs after other parallelly executed statements.
add_function(term_{1} = term_{1},
...,
term_{n} =
term_{n}).
adds reactions to compute the result of a function of the current variables in the concentration of a new variable, at steady state.
This command only accepts input list with each element of this form: <lhs> = <rhs>.
Example.
biocham:
present(x,1).biocham:
present(y,3).biocham:
add_function(z=x+y).biocham:
list_ode.[0] d(z_p)/dt=x_p+y_pz_pfast*z_p*z_n [1] d(z_n)/dt=x_n+y_nz_nfast*z_p*z_n [2] d(temp0_n)/dt=0 [3] d(temp0_p)/dt=0 [4] d(temp1_n)/dt=0 [5] d(temp1_p)/dt=0 [6] d(x)/dt=0 [7] d(x_n)/dt=0 [8] d(x_p)/dt=0 [9] d(y)/dt=0 [10] d(y_n)/dt=0 [11] d(y_p)/dt=0
biocham:
list_model.MA(1) for x_p=>z_p+x_p. MA(1) for y_p=>z_p+y_p. MA(1) for z_p=>_. MA(1) for x_n=>z_n+x_n. MA(1) for y_n=>z_n+y_n. MA(1) for z_n=>_. fast*z_p*z_n for z_n+z_p=>_. present(x,1). present(y,3). present(temp0_p,0). present(temp0_n,0). present(temp1_p,0). present(temp1_n,0). present(x_p,1). present(x_n,0). present(y_p,3). present(y_n,0). present(z_p,0). present(z_n,0). parameter( fast = 1000 ).
biocham:
numerical_simulation(time:10).biocham:
plot.Options.
 zero_ultra_mm: yesno
 set the kinetics of zeroorder ultrasensitivity to MichaelisMenten kinetics
option(simplify_variable_init:no).
option(zero_ultra_mm:no).
compile_program(Program: atom).
Using atom Program as the program code.We can write the follwing program and compile it in reactions.
biocham:
compile_program('z := 2. x := z + 1. y := x  z.',zero_ultra_mm:yes).biocham:
numerical_simulation(time:60).biocham:
plot(show:{x_p,x_n,y_p,y_n,z_p,z_n}).This is equivalent to '
add_function( z = ( 2 post ), x = ( pre z + 1 post ), y = ( pre x  z )).
' and can be compared with the parallel version of the program, without pre and post conditions.biocham:
clear_model.biocham:
compile_program('z := 2, x := z + 1, y := x  z.',zero_ultra_mm:yes).biocham:
numerical_simulation(time:20).biocham:
plot(show:{x_p,x_n,y_p,y_n,z_p,z_n}).Options.
compile_program_file(input_file).
Compile the highlevel language program from file.
Options.
11.2. Synthesis from GPAC circuits
Biocham can compile a GPAC circuit (Shannon's General Purpose Analog Computer) into an abstract CRN using positive as well as negative concentration values. Only weak GPACs, in which the integration gate is with respect to time and not a variable, are considered. The variables associated to the different points of the circuits can be named with the "::" operator. Dy default they are named x0, x1,... The syntax of weak GPAC circuits is as follows:
The option fast rate (defaulting to 1000) is intended to define a high rate constant for reactions faster than the other reactions.
It is used
 in the reactions for computing the results of the sum and product GPAC blocks,
 and in the annihilation reactions between the molecular species for the positive and negative values of a real valued variable.
option(fast_rate:1000).
compile_wgpac({wgpac_{1}, ...,
wgpac_{n}}).
compiles a set of weak GPAC circuits into a reaction network.
Options.
 fast_rate: arithmetic_expression
Example.
Compilation of a GPAC circuit generating cosine as a function of time (using positive and negative concentrations) :
biocham:
compile_wgpac(f::integral integral1*f).biocham:
present(f,1).biocham:
list_model.fast*[x2]*[f] for x2+f=>x1+x2+f. fast*[x1] for x1=>_. MA(1) for x1=>x0+x1. MA(1) for x0=>f+x0. present(x2,1). present(f,1). parameter( fast = 1000 ).
biocham:
numerical_simulation(time:10).biocham:
plot.11.3. Synthesis from polynomial differential equations
The Turing completeness of continuous CRNs [FLBP17cmsb] states that any computable function over the reals can be computed by a CRN over a finite set of molecular species. Biocham uses the proof of that result to compile any computable real function presented as the solution of a polynomial differential equation system with polynomial initial values (PIVP) into a finite CRN.
The compilation from mathematical expression is restricted to some standard functions and simple operations using a functional notation where
id
represents the operand.The compilation from PIVPs is implemented in full generality
The option for binomial reduction restricts the synthesis of reactions with at most two reactants (the default is not).
option(binomial_reduction:no).
Another option is the lazy introduction of molecular species for negative real values (the default is yes).
option(lazy_negatives:yes).
compile_from_expression(Expr: arithmetic_expression, Output: name).
creates a biochemical reaction network such that Output = Expr(time).
The expression is restricted to standard functions and simple operations using a functional notation where
id
represents the time operand.
Options.
compile_from_expression(Expr: arithmetic_expression, Input: name, Output: name).
creates a biochemical reaction network to compute the function Output = Expr(Input).
The expression is restricted to standard functions and simple operations using a functional notation where
id
represents the operand.
The Input species is initialized with the value of a special parameter named input. This species may be degraded by the computation.
Options.
Example.
Compilation of the expression 4+time^2:
biocham:
compile_from_expression(4+id*id,output).biocham:
list_model.MA(1) for A=>output+A. MA(1) for B=>output+B. MA(1) for A=>C+A. MA(1) for B=>C+B. MA(1) for _=>B. MA(1) for _=>A. present(output,4). present(D,4).
biocham:
numerical_simulation(time:4).biocham:
plot.Example.
Compilation of the expression cos(time):
biocham:
compile_from_expression(cos,costime).biocham:
list_model.fast*costime_p*costime_m for costime_p+costime_m=>_. fast*E_p*E_m for E_p+E_m=>_. MA(1) for E_p=>costime_p+E_p. MA(1) for E_m=>costime_m+E_m. MA(1) for costime_m=>E_p+costime_m. MA(1) for costime_p=>E_m+costime_p. present(costime_p,1). parameter( fast = 1000 ).
biocham:
numerical_simulation(time:10).biocham:
plot.Example.
Compilation of the expression cos(time):
biocham:
compile_from_expression(cos+1,costime1).biocham:
list_model.fast*costime1_p*costime1_m for costime1_p+costime1_m=>_. fast*C_p*C_m for C_p+C_m=>_. fast*B_p*B_m for B_p+B_m=>_. MA(1) for B_p=>costime1_p+B_p. MA(1) for B_m=>costime1_m+B_m. MA(1) for B_p=>C_p+B_p. MA(1) for B_m=>C_m+B_m. MA(1) for C_m=>B_p+C_m. MA(1) for C_p=>B_m+C_p. present(costime1_p,2). present(C_p,1). present(A,1). parameter( fast = 1000 ).
biocham:
numerical_simulation(time:10).biocham:
plot.Example.
Compilation of the expression cos(x) and computation of cos(4):
biocham:
compile_from_expression(cos,x,cosx).biocham:
list_model.fast*cosx_p*cosx_m for cosx_p+cosx_m=>_. fast*B_p*B_m for B_p+B_m=>_. MA(1) for B_p+x=>cosx_p+B_p+x. MA(1) for B_m+x=>cosx_m+B_m+x. MA(1) for cosx_m+x=>B_p+cosx_m+x. MA(1) for cosx_p+x=>B_m+cosx_p+x. MA(1) for x=>_. present(cosx_p,1). present(x,input). parameter( input = 1.0, fast = 1000 ).
biocham:
parameter(input=4).biocham:
numerical_simulation(time:10).biocham:
plot.One can also compile real valued functions defined as solutions of Polynomial Initial Value Problems (PIVP), i.e. solutions of polynomial differential equations, using the following syntax:
pode ::= 
polynomial ::= 
compile_from_pivp(PIVP: pivp, Output: name).
compiles a PIVP into a CRN (with initial concentration values) that computes the projection on one variable Output as a function of time, of the multivariate function f(t) solution of the PIVP.
Options.
 binomial_reduction: yesno
 Determine if the binomial reduction has to be performed
compile_from_pivp(PIVP: pivp, Input: name, Output: name).
compiles a PIVP into a CRN (with initial concentration values) that computes the value of variable Output at time Input, of the multivariate function f(t) solution of the PIVP.
Options.
 binomial_reduction: yesno
 Determine if the binomial reduction has to be performed
Example.
Compilation of a simple oscillator with 2 species (LotkaVolterra)
biocham:
option(lazy_negatives:yes).biocham:
compile_from_pivp((0.5,d(x)/dt=xx*y;0.25,d(y)/dt=x*y0.5*y),x).biocham:
list_model.MA(1) for x+y=>y. MA(1) for x=>2*x. y/2 for y=>_. MA(1) for x+y=>2*y+x. present(x,0.5). present(y,0.25).
biocham:
numerical_simulation(time:10).biocham:
plot(show:{x,y}).Example.
Compilation of a Hill function of order 2 as a function of input
biocham:
option(lazy_negatives:yes).biocham:
option(binomial_reduction:yes).biocham:
compile_from_pivp((0.0,d(h)/dt=2*x*x*y;1.0,d(x)/dt= 2*x*x*y;0.0,d(y)/dt=1),in,h).biocham:
parameter(input=2).biocham:
list_model.2*x*xyin for x+xyin=>h+x+xyin. MA(1) for (in)=>_. 2*x*xyin for x+xyin=>xyin. MA(1) for xin=>_. 2*xin*xyin for xin+xyin=>xyin. MA(1) for (in)+xin=>xyin+(in)+xin. MA(1) for xyin=>_. 2*xyin^2 for 2*xyin=>xyin. present(in,input). present(x,1.0). present(xin,input). parameter( input = 2 ).
biocham:
numerical_simulation(time:10).biocham:
plot(show:{h}).biocham:
plot(show:{h},logscale:x).add_reactions_from_pivp(PIVP: pivp).
creates a reaction model to implement a given PIVP (same syntax as compile_from_pivp),
variables that may became negative are splitted x => x_p  x_m, the resulting ODE
is translated in CRN through add_reactions_from_ode_system (see above).
Example.
In this example dual species are introduced for x because of its initial negative value but not for y:
biocham:
add_reactions_from_pivp((1,d(x)/dt=yx;0,d(y)/dt=2y)).biocham:
list_model.MA(1) for x_m=>x_p+x_m. MA(1) for x_p=>x_m+x_p. fast*x_p*x_m for x_m+x_p=>_. 2 for _=>y. MA(1) for y=>x_p. present(x_p,0). present(x_m,1). present(y,0). parameter( fast = 1000 ).
biocham:
numerical_simulation(time:10).biocham:
plot.11.4. Synthesis from transfer functions
transfer_polynomial ::= 
 s 
transfer_function ::= 
enable_p_m_mode.
Each variable corresponds to two species: one for the negative part and one for the positive part.
disable_p_m_mode.
Each variable corresponds to one species.
which_p_m_mode.
Displays which mode is being used.
set_p_m_rate(Rate: number).
Set the annihilation rate between two complementary species.
set_kinetics_rate(Rate: number).
Set the final summator rate
compile_transfer_function(F, U: object, Y: object).
compile a transfer function in variable s into a chemical filter.
compile_transfer_function(NL, DL, U: object, Y: object).
compile a transfer function into a chemical filter.
Example.
biocham:
set_kinetics_rate(1000).biocham:
compile_transfer_function([22,51,47,19,3],[30,66,67,36,10,1],a,b).biocham:
present(a).biocham:
numerical_simulation(method:msbdf).biocham:
plot.Chapter 12
Index
A
 about/0 1.3.
 absent/1 2.2.
 add/1 6.1.
 add_biocham/1 6.1.
 add_conservation/1 7.2.
 add_ctl/1 8.2.
 add_edge/* 3.3.
 add_event/* 5.
 add_function/* 11.1.
 add_ginml/1 6.2.
 add_influence/1 4.1.
 add_influences_from_ode/1 6.3.
 add_influences_from_ode_system/0 6.3.
 add_ode/* 6.3.
 add_qual/1 6.2.
 add_reaction/1 3.1.
 add_reactions_from_ode/1 6.3.
 add_reactions_from_ode_system/0 6.3.
 add_reactions_from_pivp/1 11.3.
 add_sbml/1 6.2.
 add_time_event/* 5.
 add_vertex/* 3.3.
 against (option) 9.3.
 alias/1 2.5.
 all_reductions (option) 7.3.
B
 bifurcations/3 9.4.
 binomial_reduction (option) 11.3.11.3.11.3.11.3.
 boolean_semantics (option) 8.1.8.2.8.2.
 boolean_simulation (option) 8.5.8.5.
C
 canonical/1 2.5.
 change_parameter_to_variable/1 9.4.
 check_conservations/0 7.2.
 check_ctl/0 8.2.
 check_multistability/0 7.6.
 check_oscillations/0 7.6.
 cleanup_ctl/0 8.2.
 clear_dimension/1 7.1.
 clear_dimensions/0 7.1.
 clear_initial_state/0 2.2.
 clear_model/0 6.1.
 clear_variation/1 9.4.
 clear_variations/0 9.4.
 cmaes_init_center (option) 10.3.10.3.
 cmaes_log_normal (option) 10.3.10.3.
 cmaes_stop_fitness (option) 10.3.10.3.
 cnf_clause_size (option) 8.5.8.5.
 compile_from_expression/2 11.3.
 compile_from_expression/3 11.3.
 compile_from_pivp/2 11.3.
 compile_from_pivp/3 11.3.
 compile_program/1 11.1.
 compile_program_file/1 11.1.
 compile_transfer_function/3 11.4.
 compile_transfer_function/4 11.4.
 compile_wgpac/1 11.2.
 conservation_size (option) 7.2.
 continue/1 9.1.
D
 data_transform (option) 8.5.8.5.
 delete/* 6.1.
 delete_alias/1 2.5.
 delete_attribute/2 3.3.
 delete_column/* 10.1.
 delete_conservation/1 7.2.
 delete_conservations/0 7.2.
 delete_ctl/0 8.2.
 delete_ctl/1 8.2.
 delete_edge/* 3.3.
 delete_function/* 2.4.
 delete_graph/1 3.3.
 delete_ltl_pattern/* 10.2.
 delete_ode/* 6.3.
 delete_ode_system/1 6.3.
 delete_parameter/* 2.3.
 delete_reaction/1 3.1.
 delete_reaction_named/1 3.1.
 delete_reaction_prefixed/1 3.1.
 delete_reactions/0 3.1.
 delete_row/* 10.1.
 delete_table/* 10.1.
 delete_temporary_files/0 1.3.
 delete_vertex/* 3.3.
 disable_p_m_mode/0 11.4.
 distinct_species (option) 7.3.
 dose_response/3 9.4.
 double_michaelis_menten (option) 7.4.
 download_curated_biomodel/1 6.2.
 draw_first (option) 3.2.
 draw_graph/0 3.3.
 draw_influence_hypergraph/0 4.2.
 draw_influences/0 4.2.
 draw_neighborhood/0 2.2.
 draw_reactions/0 3.2.
 draw_rules/0 3.2.
E
 enable_p_m_mode/0 11.4.
 enzyme (option) 7.4.
 ep (option) 7.4.
 error_epsilon_absolute (option) 9.1.10.3.10.3.10.3.
 error_epsilon_relative (option) 9.1.10.3.10.3.10.3.
 expand_ctl/1 8.2.
 expand_ltl/1 10.2.
 export_biocham/1 6.1.
 export_graph/1 3.3.
 export_lemon_graph/1 4.2.
 export_nusmv/1 8.2.
 export_ode/1 6.3.
 export_ode_as_latex/1 6.3.
 export_plot/1 9.3.
 export_plot_to_canvas/2 9.3.
 export_plot_to_png/1 9.3.
 export_sbml/1 6.2.
 export_table/1 10.1.
 extremal_sepi (option) 7.3.
F
 fast_rate (option) 11.2.
 filter (option) 9.1.
 find_reaction_prefixed/2 3.1.
 foltl_magnitude (option) 10.2.10.2.
 force_graph (option) 4.2.
 function/* 2.4.
G
H
 hill_reaction (option) 7.4.
 hybrid_dynamic_simulation/3 9.2.
 hybrid_dynamic_simulation/6 9.2.
 hybrid_static_simulation/2 9.2.
 hybrid_static_simulation/5 9.2.
I
 import_ode/1 6.3.
 import_reactions_from_graph/0 3.2.
 import_reactions_with_inhibitors (option) 6.3.6.3.6.3.6.3.
 influence_graph/0 4.2.
 influence_hypergraph/0 4.2.
 influence_model/0 4.1.
 inherits/* 6.1.
 init/1 6.3.
 initial_step_size (option) 9.1.10.3.10.3.10.3.
K
 keep_temporary_files/0 1.3.
L
 lazy_negatives (option) 11.3.11.3.
 left_to_right (option) 3.3.
 list_aliases/0 2.5.
 list_attributes/1 3.3.
 list_columns/0 10.1.
 list_conservations/0 7.2.
 list_ctl/0 8.2.
 list_current_models/0 6.1.
 list_dimensions/0 7.1.
 list_edges/0 3.3.
 list_events/0 5.
 list_functions/0 2.4.
 list_graph_objects/0 3.3.
 list_graphs/0 3.3.
 list_influences/0 4.1.
 list_initial_state/0 2.2.
 list_isolated_vertices/0 3.3.
 list_last_state/0 10.1.
 list_locations/0 2.2.
 list_ltl_patterns/0 10.2.
 list_model/0 6.1.
 list_models/0 6.1.
 list_molecules/0 2.2.
 list_neighborhood/0 2.2.
 list_ode/0 6.3.
 list_ode_systems/0 6.3.
 list_options/0 1.3.
 list_parameters/0 2.3.
 list_reactions/0 3.1.
 list_rows/0 10.1.
 list_rules/0 3.1.
 list_stable_states/0 8.1.
 list_tables/0 10.1.
 list_tscc_candidates/0 8.1.
 load/1 6.1.
 load_biocham/1 6.1.
 load_ginml/1 6.2.
 load_influences_from_ode/1 6.3.
 load_influences_from_ode_system/0 6.3.
 load_qual/1 6.2.
 load_reactions_from_ode/1 6.3.
 load_reactions_from_ode_system/0 6.3.
 load_sbml/1 6.2.
 load_table/1 10.1.
 logarithmic_bifurcations/3 9.4.
 logarithmic_dose_response/3 9.4.
 logarithmic_variation/3 9.4.
 logscale (option) 9.3.
 ltl_pattern/* 10.2.
M
 make_absent_not_present/0 2.2.
 make_present_not_absent/0 2.2.
 mapping_restriction (option) 7.3.
 max_nb_reductions (option) 7.3.
 maximum_step_size (option) 9.1.10.3.10.3.10.3.
 merge_reactions/2 3.1.
 merge_restriction (option) 7.3.
 method (option) 9.1.10.3.10.3.10.3.
 michaelis_menten (option) 7.4.
 michaelis_menten_expansion (option) 7.4.
 minimum_step_size (option) 9.1.
 multistability_graph/0 4.2.
N
 new_graph/0 3.3.
 new_model/0 6.1.
 new_ode_system/0 6.3.
 numerical_simulation/0 9.1.
 nusmv_counter_example (option) 8.2.
 nusmv_initial_states (option) 8.2.
O
P
 pac_learning/1 8.5.
 pac_learning/3 8.5.
 parameter/* 2.3.
 parametrize/0 6.1.
 partial_hill_reaction (option) 7.4.
 pattern_reduction/2 7.4.
 place/* 3.3.
 plot/0 9.3.
 precision (option) 9.1.10.3.10.3.10.3.
 present/1 2.2.
 present/2 2.2.
 prolog/1 1.3.
Q
R
 r_1 (option) 7.4.
 r_2 (option) 7.4.
 reaction_graph/0 3.2.
 reaction_model/0 4.1.
 reduce_model/0 8.3.
 reduce_model/1 8.3.
 rename_column/2 10.1.
 rename_table/1 10.1.
 reset_options/0 1.3.
 restore_parameter/1 9.4.
 revise_model/0 8.4.
 revise_model/1 8.4.
 robustness/3 10.3.
 robustness_coeff_var (option) 10.3.
 robustness_relative_error (option) 10.3.
 robustness_samples (option) 10.3.
 rule_graph/0 3.2.
S
 satisfaction_degree/2 10.2.
 search_conservations/0 7.2.
 search_parameters/2 10.3.
 search_parameters/3 10.3.
 search_reduction/2 7.3.
 seed/1 1.3.
 select_graph/1 3.3.
 select_model/1 6.1.
 select_ode_system/1 6.3.
 select_table/1 10.1.
 set_attribute/2 3.3.
 set_dimension/3 7.1.
 set_graph_name/1 3.3.
 set_kinetics_rate/1 11.4.
 set_model_name/1 6.1.
 set_ode_system_name/1 6.3.
 set_p_m_rate/1 11.4.
 show (option) 9.3.9.4.9.4.9.4.9.4.
 simplify_variable_init (option) 11.1.11.1.
 stats (option) 7.3.7.5.9.1.
 stochastic_conversion (option) 9.1.
 stochastic_thresholding (option) 9.1.
T
 test_permutations (option) 7.6.
 test_transpositions (option) 7.6.
 time (option) 9.1.9.4.9.4.9.4.9.4.9.4.9.4.10.3.10.3.10.3.
 timeout (option) 7.3.
 transition/* 3.3.
 tropical_denominator (option) 7.5.
 tropical_epsilon (option) 7.5.
 tropical_ignore (option) 7.5.
 tropical_max_degree (option) 7.5.
 tropical_single_solution (option) 7.5.
 tropicalize/0 7.5.
U
 undefined/1 2.2.
V
W
 which_p_m_mode/0 11.4.