logo.pngBiocham 4.4.4 Reference Manual

Inria http://lifeware.inria.fr/biocham4/

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 :
Biocham v4 is a complete rewriting (in SWI-Prolog) of Biocham v3 (written in GNU-Prolog). It introduces some new features, mainly:
plus since v4.1:
Some features of Biocham v3 are still not implemented in v4, mainly:
or less efficiently (because of on-the-fly C compilation instead of native Prolog code):

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:
Biocham numerical data time series can be imported⁄exported in
The following files can also be used to export some Biocham objects:

Biocham call options

biocham [file] launches Biocham and optionnally loads the file given as argument.
biocham –trace –version –list-commands 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://jupyter-notebook.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 notebook-kernel communication timeout of 120 (seconds).
To execute a Biocham command enter Shift-Return. 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 Biocham-4 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 top-level

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.0e-6)
[38] option(error_epsilon_relative:1.0e-6)
[39] option(initial_step_size:1.0e-6)
[40] option(maximum_step_size:0.01)
[41] option(minimum_step_size:1.0e-5)
[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=b-a
[1] d(b)/dt=a-b
biocham: option(time:5).

biocham: numerical_simulation(method:ssa). plot.
Simulation results
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:
Both types of models can be combined and interpreted in a hierarchy of semantics, including the differential semantics (Ordinary Differential Equations), stochastic semantics (Continuous-time 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 ::=
| atom
output_file ::=
| atom
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 ::=
| atom
identifier ::=
| name @ name
| name
parameter_name ::=
function_prototype ::=
| term
object ::=
arithmetic_expression ::=
| [ object ]
| product( term in term , arithmetic_expression )
condition ::=
| true
| false
| not condition
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({object1, ..., objectn}).
Every object in [object1, ..., objectn] is made present in the initial state with initial concentration equal to 1.
present({object1, ..., objectn}, concentration).
Every object in [object1, ..., objectn] 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({object1, ..., objectn}).
Every object in [object1, ..., objectn] is made absent in the initial state.
undefined({object1, ..., objectn}).
Every object in [object1, ..., objectn] 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_name1 = number1, ..., parameter_namen = numbern).
sets the value of parameters.
list_parameters.
shows the values of all known parameters.
delete_parameter(parameter_name1, ..., parameter_namen).
deletes some parameters

2.4. Functions

function(function_prototype1 = arithmetic_expression1, ..., function_prototypen = arithmetic_expressionn).
sets the definition of functions.
list_functions.
lists all known functions.
delete_function(functor1, ..., functorn).
deletes some functions. Either arity is given, or all functions with the given functor are deleted.

2.5. Aliases

alias(object1 = ... = objectn).
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 ::=
| name
basic_reaction ::=
reactants ::=
solution ::=
| _
kinetics ::=
Useful abbreviations for mass action law kinetics (with inhibitors), Michaelis-Menten 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(reaction1, reaction2).
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.
Simulation results
reaction_graph.
Builds the reaction graph of the current model.

Example.
biocham: option(draw_first:none,left_to_right:no).

biocham: draw_reactions.
Simulation results
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

attribute ::=
| name = term
| name
edge ::=
| name -> name
edgeref ::=
| edge
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(name1, ..., namen).
Adds new vertices to the current graph.
delete_vertex(name1, ..., namen).
Deletes a set of vertices from the current graph.
add_edge(edge1, ..., edgen).
Adds the given set of edges to the current graph. The vertices are added if needed.
delete_edge(edgeref1, ..., edgerefn).
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.
graph_object ::=
| edge
| name
set_attribute({graph_object1, ..., graph_objectn}, attribute).
Adds an attribute to every vertex or edge in the given set. The vertices and the edges are added if needed.
place(name1, ..., namen).
Sets that the vertices [name1, ..., namen] are places.
transition(name1, ..., namen).
Sets that the vertices [name1, ..., namen] are transitions.
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 ::=
basic_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.
Simulation results
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 molecule-influence 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.
Simulation results
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 differential-stochastic simulators [CFJS15tomacs].
add_event(condition, parameter_name1 = arithmetic_expression1, ..., parameter_namen = arithmetic_expressionn).
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.
Simulation results
biocham: numerical_simulation(time:2,method:ssa).

biocham: plot.
Simulation results
Note that the continuous numerical_simulation/0 engine will attempt to interpolate linear event conditions as per [doi.org/10.1007/3-540-45351-2_19].
add_time_event(condition, parameter_name1 = arithmetic_expression1, ..., parameter_namen = arithmetic_expressionn).
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

range ::=
ref ::=
| [ range , ... , range ]
| name
load(input_file).
acts as the corresponding load_biocham/1load_sbml/1load_ode/1load_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/1add_sbml/1add_ode/1add_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({ref1, ..., refn}).
selects some models.
set_model_name(name).
changes the current model name.
delete([range1], ..., [rangen]).
deletes the listed elements from the model.
inherits(ref1, ..., refn).
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 SBML-qual 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 ::=
| name

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(ode1, ..., oden).
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(oderef1, ..., oderefn).
removes the given set of ODEs from the current ODE system.
init(name1 = arithmetic_expression1, ..., namen = arithmetic_expressionn).
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 Place-invariants 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 P-invariants of the network up to the maximal size given by the option conservation_size and defaulting to 4. Such P-invariants are particular mass conservation laws that are independent from the kinetics.
Options.
conservation_size: integer
Set the maximal stoichiometric size for a P-invariant.

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 NP-complete problem.
extremal_sepi ::=
| no
Search standard reduction.
| minimal_deletion
Search reduction that minimises bottom (i.e., the number of deletions).
| maximal_deletion
Search reduction that maximises bottom (i.e., the number of deletions).
merge_restriction ::=
| no
Search standard reduction.
| neighbours
Search reduction with local two-neighbour restriction (cf. report chapter 5).
| old
Search reduction with old two-neighbour restriction (prior to 2019).
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_influence1, ..., basic_influencen]
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 Michaelis-Menten 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 Michaelis-Menten 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 Michaelis-Menten patterns, in order to exhibit the intermediary step of reversible complexation enzyme-substrate.
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 ::=
| yes
Value yes.
| no
Value no.
| maybe
Value 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 one-reaction Michaelis-Menten pattern, with a double-arrow from the enzyme. The michaelis_menten_expansion option has the opposite effect to expand reduced Michaelis-Menten patterns
Options.
michaelis_menten: yesno
specifies if reducing Michaelis-Menten patterns
r_1: yesnomaybe
R-1 reaction (from the commplex ES to enzyme E and substrate S) in the searched Michaelis-Menten pattern
r_2: yesnomaybe
R-2 reaction (from the product P and enzyme E to the complex EP or ES) in the searched Michaelis-Menten pattern
ep: yesnomaybe
irreversible reaction from the complex ES to the complex EP in the searched Michaelis-Menten 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 Michaelis-Menten
enzyme: yesno
specifies if the reduced Michaelis-Menten pattern (that will appear in the output graph) is with enzyme or no
michaelis_menten_expansion: yesno
expansion of Michaelis-Menten patterns (form the reduced form with enzyme to form with complex ES and reaction R-1), should be used with option michaelis_menten:no to avoid confusion

Example. Example of reduction of a Michaelis-Menten pattern with reactions R-1, R-2 and complex EP. Input graph :
reaction graph with the reversible formation of the complex enzyme-substrate ES, the irreversible transformation of the complex from ES to EP and the reversible decomplexation from EP to enzyme and product
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 Michaelis-Menten pattern) :
Reduced Michaelis-Menten reaction

Example. Hill pattern (can be reduced with option hill_reaction) :
Hill reaction (two ligant sites on the enzyme)

Example. Partial Hill pattern (can be reduced with option partial_hill_reaction) :
Partial Hill reaction (twice the same ligant site on the enzyme)

Example. Double Michaelis-Menten pattern (can be reduced with option double_michaelis_menten) :
Double Michaelis-Menten reaction (two forms of the complex, ES and ESS, can release the product)

Example. Example of expansion of reduced Michaelis-Menten in a MAPK cascade.
biocham: load(library:examples/sepi/mapk0.bc).

biocham: draw_reactions.
Simulation results
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.
Simulation results

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 fast-slow dynamics leading to model reductions based on quasi-steady states or reactions in quasi-equilibrium [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+Cdc2-Cyclin~{p1,p2}+Cdc2-Cyclin~{p1}+Cdc2~{p1}
1 complex invariant(s)

found a complete equilibration leading to the rescaling:
Cdc2' = ε^(- 2) * Cdc2
Cdc2-Cyclin~{p1,p2}' = Cdc2-Cyclin~{p1,p2}
Cdc2-Cyclin~{p1}' = ε^(- 2) * Cdc2-Cyclin~{p1}
Cdc2~{p1}' = Cdc2~{p1}
Cyclin' = ε^(- 4) * Cyclin
Cyclin~{p1}' = ε^(- 2) * Cyclin~{p1}


found a complete equilibration leading to the rescaling:
Cdc2' = ε^(- 3) * Cdc2
Cdc2-Cyclin~{p1,p2}' = Cdc2-Cyclin~{p1,p2}
Cdc2-Cyclin~{p1}' = ε^(- 2) * Cdc2-Cyclin~{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
Cdc2-Cyclin~{p1,p2}' = Cdc2-Cyclin~{p1,p2}
Cdc2-Cyclin~{p1}' = ε^(- 2) * Cdc2-Cyclin~{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: {object1, ..., objectn}
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 non-degenerate 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 non-degenerate (i.e. with no variable equal to 0) steady states.
Options.
test_transpositions: yesno
Lazy test of all possible swappings between two species only
test_permutations: yesno
Lazy test of all possible permutations between species. This option may highly increase the computation time.

Example.
biocham: load("library:biomodels/BIOMD0000000040.xml").

biocham: multistability_graph.

biocham: draw_graph(left_to_right:yes).
Simulation results
biocham: check_multistability.
There may be non-degenerate multistationarity, positive circuit detected.
biocham: check_multistability(test_transpositions:yes).
There is no multiple steady states with non-zero 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 ::=
| EX( ctl )
| EF( ctl )
| EG( ctl )
| EU( ctl , ctl )
| ER( ctl , ctl )
| AX( ctl )
| AF( ctl )
| AG( ctl )
| AU( ctl , ctl )
| AR( ctl , ctl )
| not ctl
| ctl /\ ctl
| ctl \/ ctl
| ctl -> ctl
| reachable( 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 model-checkers. 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 some
option(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 model-checker. As is usual in Model-Checking, 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 counter-example 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 'RAF-RAFK'=>RAF+RAFK,MA(1) for 'RAF~{p1}'+RAFPH=>'RAF~{p1}-RAFPH',MA(1) for 'RAF~{p1}-RAFPH'=>'RAF~{p1}'+RAFPH,MA(1) for 'MEK-RAF~{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 'MAPK-MEK~{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',1e-5,1e-3, time:200, show:'MAPK~{p1,p2}').
Simulation results

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 theory-revision on the current model to satisfy the query given as argument.
revise_model.
Use theory-revision as above, using the currently defined CTL specification.

8.5. PAC Learning influence models from Boolean traces

Implementation of Leslie Valiant's PAC-Learning [doi.org/10.1215/0961754X-2872666] algoritm for k-CNF 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
| threshold( integer )
generate_traces(NInitialStates: integer, timeHorizon, FilePrefix: output_file).
Randomly choses NInitialStates initial states and for each, runs a numerical simulation with method: spn, method: ssa or method: sbn and timeHorizon 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.
data_transform: transform
transformation of the traces before saving
boolean_simulation: yesno
use a boolean simulation (sbn) instead of the default stochastic one
pac_learning({input_file1, ..., input_filen}).
Uses every time-series .csv file in [input_file1, ..., input_filen] 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 Lotka-Volterra 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.
cnf_clause_size: integer
Maximum size of CNF clauses learnt.
data_transform: transform
transformation of the traces before learning
boolean_simulation: yesno
use a boolean simulation (sbn) instead of the default stochastic one (spn) to generate the data


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/Stepping-Functions.html#Stepping-Functions 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 ::=
| rk2
Explicit embedded Runge-Kutta (2, 3) method
| rk4
Explicit 4th order (classical) Runge-Kutta. Error estimation is carried out by the step doubling method. For more efficient estimate of the error, use the embedded methods described below
| rkf45
Explicit embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator
| rkck
Explicit embedded Runge-Kutta Cash-Karp (4, 5) method.
| rk8pd
Explicit embedded Runge-Kutta Prince-Dormand (8, 9) method.
| rk1imp
Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method.
| rk2imp
Implicit Gaussian second order Runge-Kutta. Also known as implicit mid-point rule. Error estimation is carried out by the step doubling method.
| rk4imp
Implicit Gaussian 4th order Runge-Kutta. Error estimation is carried out by the step doubling method.
| bsimp
Default method. Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.
| msadams
A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12.
| msbdf
Perhaps the most robust method but may be slow, a variable-coefficient 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 non-linear equations. Method order varies dynamically between 1 and 5.
| ssa
Stochastic simulation of a Continuous-Time 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.
| spn
Stochastic Petri net simulation. Similar to the SSA algorithm above but with a discrete⁄logical time.
| pn
Random Petri net simulation run. All transitions are equiprobable.
| sbn
Stochastic Boolean net simulation. Similar to the SSA algorithm above but with a discrete⁄logical time and Boolean states.
| bn
Random Boolean net simulation run. All transitions are equiprobable.
| rsbk
Rosenbrock method.
time ::=
filter ::=
| no_filter
Does no filtering at all.
| only_extrema
Only keeps the points that are an extremum for some variable.
option(method:bsimp).
option(error_epsilon_absolute:1.0e-6).
option(error_epsilon_relative:1.0e-6).
option(initial_step_size:1.0e-6).
option(maximum_step_size:0.01).
option(minimum_step_size:1.0e-5).
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 user-defined functions over time.

Example.
biocham: load(library:examples/mapk/mapk).

biocham: numerical_simulation(method:msbdf).

biocham: plot.
Simulation results
biocham: plot(show: {MAPK~{p1,p2}, MEK~{p1,p2}}).
Simulation results
biocham: numerical_simulation(filter:only_extrema).

biocham: plot.
Simulation results
Options.
logscale: axes
Apply log-scaling to the specified axes.
show: {object1, ..., objectn}
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.
Simulation results
biocham: numerical_simulation(time:5,method:msbdf).

biocham: plot.
Simulation results
biocham: plot(logscale:y).
Simulation results
biocham: plot(show:a,against:b).
Simulation results
biocham: plot(show:a,against:b,logscale:xy).
Simulation results
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 dose-response 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.0e-7,0.1).

biocham: logarithmic_variation(b,1.0e-7,0.1).

biocham: numerical_simulation.

biocham: plot.
Simulation results
biocham: plot(logscale:y).
Simulation results
biocham: list_model.
variation_a--0.016666649999999998*1 for _=>a.
variation_b--b/6*log(999999.9999999) for b=>2*b.
present(a,1.0e-7).
present(b,1.0e-7).
dose_response(Dose: object, Value1: arithmetic_expression, Value2: arithmetic_expression).
Draws a dose-response 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.
time: time
time horizon necessary to approximate the steady state response; the variation of the dose is performed over a time horizon of 10 fold that value.
show: {object1, ..., objectn}
Defines a set of objects for the response; everything will be plotted if the show option is empty
logarithmic_dose_response(Dose: object, Value1: arithmetic_expression, Value2: arithmetic_expression).
Draws a similar dose-response diagram by logarithmic variation and plots it on a logarithmic scale for both the dose and the responses. Value1 must be non zero.
Options.
time: time
show: {object1, ..., objectn}

Example. This example shows the dose-response diagrams of the MAPK model of Huang and Ferrell in BiomModels revealing the amplifier and analog-digital 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,1e-6,1e-4).
Simulation results
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 dose-response 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.
Options.
time: time
show: {object1, ..., objectn}
logarithmic_bifurcations(Dose: object, Value1: arithmetic_expression, Value2: arithmetic_expression).
Draws a similar dose-response diagram by logarithmic variation using a logarithmic scale for both the dose and the response. Value1 and Value2 must be non zero.
Options.
time: time
show: {object1, ..., objectn}

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).
Simulation results
biocham: bifurcations(MAPKK, 0, 100, time:1e4, show:Mpp).
Simulation results

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(name1, ..., namen).
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.
column ::=
| name
delete_column(column1, ..., columnn).
deletes the given columns from the current table.
rename_column(column, name).
renames the given column of the current table.
delete_row(number1, ..., numbern).
deletes the given rows from the current table.
list_last_state.
lists the last state of the current simulation trace.

10.2. Temporal logic FO-LTL(Rlin) formulae

Firt-Order Linear Time Logic with linear constraints over the reals, FO-LTL(Rlin), can be used to specify semi-qualitative semi-quantitative constraints on the dynamical behavior of the system, in a much more flexible manner than by curves to fit [RBFS11tcs]. The syntax of FO-LTL(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 ::=
| X( foltl )
| F( foltl )
| G( foltl )
| not foltl
| U( foltl , foltl )
| R( foltl , foltl )
| W( foltl , foltl )
| foltl /\ foltl
| foltl \/ foltl
| exists( name , foltl )
| foltl => foltl
| foltl <=> foltl
| false
| true
| name(untyped , ... , untyped)
foltl_predicate ::=
time_value ::=
time_value_list ::=
| [ time_value , ... , time_value ]
parameter_between_interval ::=
variable_objective ::=
| name -> number
Furthermore, some useful abbreviations have been predefined with the following FO-LTL(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=Maximum-Minimum)).
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=T2-T1)))).
ltl_pattern first_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, first_successive_peaks(Molecule, C1, T1, C2, T2) /\ Period=T2-T1)))).
ltl_pattern last_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, last_successive_peaks(Molecule, C1, T1, C2, T2) /\ Period=T2-T1)))).
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=T2-T1))))).
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=T2-T1)))))).
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=T2-T1)))).
ltl_pattern base_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, successive_bases(Molecule, C1, T1, C2, T2) /\ Period=T2-T1)))).
ltl_pattern first_base_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, first_successive_bases(Molecule, C1, T1, C2, T2) /\ Period=T2-T1)))).
ltl_pattern last_base_period(Molecule, Period) = exists(T1, exists(T2, exists(C1, exists(C2, last_successive_bases(Molecule, C1, T1, C2, T2) /\ Period=T2-T1)))).
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=T2-T1))))).
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=T2-T1)))))).
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=T2-T1)))).
FO-LTL(Rlin) formulas are evaluated on traces, i.e. numerical data time series either generated by simulation or loaded from biological experiments. FO-LTL(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.
Simulation results
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_objective1, ..., variable_objectiven]).
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_prototype1 = foltl1, ..., function_prototypen = foltln).
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(functor1, ..., functorn).
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. FO-LTL(Rlin) properties

The continuous satisfaction degree of an FO-LTL(Rlin) in a given trace with respect to the objective values for the free variables can be used to compute parameter sensitivity indices and robustness measures with respect to parameter perturbation according to normal distributions, and to search parameter values for satisfying an FO-LTL(Rlin) formula [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_name1, ..., parameter_namen], [variable_objective1, ..., variable_objectiven]).
computes the robustness degree as defined in [RBFS11tcs], with respect to formula Formula, for the list of parameters [parameter_name1, ..., parameter_namen] and with list of objectives for the free variables of Formula given in [variable_objective1, ..., variable_objectiven]. 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_interval1, ..., parameter_between_intervaln], [variable_objective1, ..., variable_objectiven]).
tries to satisfy a FO-LTL(Rlin) constraint by varying the parameters listed in [parameter_between_interval1, ..., parameter_between_intervaln]. This search command uses the stochastic optimization procedure CMA-ES (covariance matrix adaptation evolution strategy) with the continuous satisfaction degree (truncated to 1 or not according to stop option) of the given FO-LTL(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.74e-01 <= stopFitness (1.00e-04)
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.98e-01 <= stopFitness (1.00e-04)
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_interval1, ..., parameter_between_intervaln], 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(a-b<x/\b-a<x), [], [x -> 10]), (G(a-b<y/\b-a<y), [kb=2], [y -> 10])]).
Time: 5.459 s
Stopping reason: Fitness: function value 4.40e-05 <= stopFitness (1.00e-04)
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 ::=
| name
expression ::= The support of arithmetic_expression is not complete, but most of them are supported.
control_flow_statement ::=
| if condition
| while( condition )
| else
| endif
| endwhile
Statements terminated by ',' will execute parallelly with next statement.
lhs ::=
rhs ::=
| para expression
| pre expression
| pre expression post
| expression post
control_flow_keyword ::=
| if_tag
| while
When using these two keywords, the rhs serves as the condition.
| else_tag
| endif
| endwhile
When using these three keywords, the rhs does not matter.
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(term1 = term1, ..., termn = termn).
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_p-z_p-fast*z_p*z_n
[1] d(z_n)/dt=x_n+y_n-z_n-fast*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.
Simulation results
Options.
zero_ultra_mm: yesno
set the kinetics of zero-order ultrasensitivity to Michaelis-Menten 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}).
Simulation results
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}).
Simulation results
Options.
simplify_variable_init: yesno
set if the variable initialization in program need to be simplified to present command
zero_ultra_mm: yesno
set the kinetics of zero-order ultrasensitivity to Michaelis-Menten kinetics
compile_program_file(input_file).
Compile the high-level language program from file.
Options.
simplify_variable_init: yesno
set if the variable initialization in program need to be simplified to present command
zero_ultra_mm: yesno
set the kinetics of zero-order ultrasensitivity to Michaelis-Menten kinetics

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:
wgpac ::=
wgpac_box ::=
| integral wgpac
wgpac_named ::=
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
option(fast_rate:1000).
compile_wgpac({wgpac1, ..., wgpacn}).
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 integral-1*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.
Simulation results

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.
binomial_reduction: yesno
Determine if the binomial reduction for synthesizing reactions with at most two reactants has to be performed
lazy_negatives: yesno
Switch between systematic and lazy introduction of molecular species for negative values
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.
binomial_reduction: yesno
Determine if the binomial reduction for synthesizing reactions with at most two reactants has to be performed
lazy_negatives: yesno
Switch between systematic and lazy introduction of molecular species for negative values

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.
Simulation results

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.
Simulation results

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.
Simulation results

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.
Simulation results
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:
pivp ::=
| pivp ; pivp
pode ::=
| d( name )/dt= polynomial
polynomial ::=
monomial ::=
| name
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 (Lotka-Volterra)
biocham: option(lazy_negatives:yes).

biocham: compile_from_pivp((0.5,d(x)/dt=x-x*y;0.25,d(y)/dt=x*y-0.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}).
Simulation results

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}).
Simulation results
biocham: plot(show:{h},logscale:x).
Simulation results
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=y-x;0,d(y)/dt=2-y)).

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.
Simulation results

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.
Simulation results

Chapter 12
Index

A

B

C

D

E

F

G

H

I

K

L

M

N

O

P

Q

R

S

T

U

V

W

X

Y

Z