logo.pngBiocham 4.6.26 Reference Manual

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

Chapter 1
Overview

1.1. Main features

help(Command: name).
Provides online help about its Command argument.
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, call options and notebook

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 —-version —-list-commands —-trace are possible options (prefixed by two dashes).
biocham_debug is the command to use to launch Biocham with the Prolog toplevel.

Jupyter notebook with Biocham kernel

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 new biocham notebook can be created with the Jupyter menu "new". The examples section of Biocham server and directory contains several notebooks.
To execute a Biocham command in a Jupyter notebook enter Shift-Return. All shortcuts are described in the keyboard menu and automatic completion is active.
The only magic commands available in the notebook and not in Biocham are the following three directives:
%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.
%timeout s, which modifies the BIOCHAM_TIMEOUT environment variable to the prescribed value (in sec.). This command modifies the timeout variable for the whole notebook and should be called in an independent cell.
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).
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 former Biocham GUI based on Jupyter is no longer available.

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) can take named options as arguments.
All the options can be defined either globally for the whole model with the command option(Option: Value, ..., Option: Value), or locally for a single command by adding additional arguments of the form Option: Value. Local options take precedence over global options.
Caution, the command clear_model does not change the current options and there is no command to restore the default options.
list_options.
lists the set of options defined in the current model with their current values (default values below).

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(use_sbml_id:no)
[4] option(output_to_library:no)
[5] option(import_reactions_with_inhibitors:yes)
[6] option(second_order_closure:no)
[7] option(infer_hidden_molecules:no)
[8] option(conservation_size:4)
[9] option(mapping_restriction:[])
[10] option(merge_restriction:no)
[11] option(timeout:180)
[12] option(all_reductions:no)
[13] option(distinct_species:no)
[14] option(max_nb_reductions:200)
[15] option(extremal_sepi:no)
[16] option(show_support:no)
[17] option(michaelis_menten:yes)
[18] option(r_1:yes)
[19] option(r_2:no)
[20] option(ep:no)
[21] option(enzyme:yes)
[22] option(hill_reaction:no)
[23] option(partial_hill_reaction:no)
[24] option(double_michaelis_menten:no)
[25] option(michaelis_menten_expansion:no)
[26] option(tropical_epsilon:0.1)
[27] option(tropical_max_degree:3)
[28] option(tropical_ignore:[{}])
[29] option(tropical_denominator:0)
[30] option(tropical_single_solution:no)
[31] option(test_permutations:no)
[32] option(test_transpositions:no)
[33] option(boolean_semantics:negative)
[34] option(boolean_initial_states:all)
[35] option(boolean_trace:no)
[36] option(nusmv_topological_order:no)
[37] option(query:current_spec)
[38] option(boolean_state_display:present)
[39] option(data_transform:none)
[40] option(cnf_clause_size:3)
[41] option(boolean_simulation:no)
[42] option(method:bsimp)
[43] option(error_epsilon_absolute:1.0e-6)
[44] option(error_epsilon_relative:1.0e-6)
[45] option(initial_step_size:1.0e-6)
[46] option(maximum_step_size:0.01)
[47] option(minimum_step_size:1.0e-5)
[48] option(precision:6)
[49] option(time:20)
[50] option(steps:1)
[51] option(stochastic_conversion:100)
[52] option(stochastic_bound:1000000.0)
[53] option(stochastic_thresholding:1000)
[54] option(filter:no_filter)
[55] option(stats:no)
[56] option(plot_table:'')
[57] option(logscale:'')
[58] option(against:Time)
[59] option(xmin:auto)
[60] option(ymin:auto)
[61] option(xmax:auto)
[62] option(ymax:auto)
[63] option(foltl_precision:12)
[64] option(foltl_magnitude:5)
[65] option(robustness_coeff_var:0.1)
[66] option(robustness_samples:30)
[67] option(positive_sampling:yes)
[68] option(openmpi_procs:1)
[69] option(cmaes_init_center:no)
[70] option(cmaes_log_normal:no)
[71] option(cmaes_stop_fitness:0.0001)
[72] option(resolution:10)
[73] option(simplify_variable_init:no)
[74] option(zero_ultra_mm:no)
[75] option(quadratic_reduction:fastnSAT)
[76] option(quadratic_bound:carothers)
[77] option(lazy_negatives:yes)
[78] option(negation_first:yes)
[79] option(fast_rate:1000)
[80] 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
[2] init(a=1)
[3] init(b=0)





biocham: option(time:5).






biocham: numerical_simulation(method:ssa). plot.


Simulation results











option.


sets options globally.




  
Example.

    


biocham: option(time:20, method:bsimp).





     









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. E.g. when you need to modify the SWI-prolog flags from
    the toplevel of Biocham: "prolog('set_prolog_flag(stack_limit, 2 000 000 000)')."



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.



with_timer(Goal, Time).


Execute a given Goal in a given Time (in ms). If no second argument is given, it displays the time taken to execute the command in the standard output.



with_timer(Goal).









  
Example.




biocham: with_timer(a+b => c).


The previous command takes 0.777006 ms











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
Basic Syntax
2.1. Names and arithmetic expressions


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


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.



2.2. Initial state


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.



initial_state(object1 = concentration1,
       ...,
       objectn =
       concentrationn).


Sets the value of initial concentration.



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 = arithmetic_expression1,
       ...,
       parameter_namen =
       arithmetic_expressionn).


sets parameter values to numbers or to the result of ground arithmetic expressions (not involving user defined functions or other 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.



evaluate(Expr: arithmetic_expression).


Evaluates an arithmetic expression possibly containing user-defined functions, parameters and molecule names denoting their initial concentration values. Warning: does loop in case of cyclic definitions.



Chapter 3
Reaction Networks
3.1. Syntax of reactions


(Bio)chemical reaction networks (CRNs) can be described in BIOCHAM by a reaction model, 
i.e. a finite multiset of directed reaction rules, formed with reactants, products, catalysts and possibly inhibitors, using the syntax defined by the grammar below.
Catalysts are molecular species which is both reactant and product in a reaction.
They can be noted with an abbreviated syntax between brackets on the arrow.
The inhibitors are distinguished in the reactants by following the symbol "⁄"  in the left-hand side of the reactions.
When not specified, the rate function is by default a mass action law kinetics with rate constant equal to 1.
BIOCHAM reaction models are compatible with the Systems Biology Markup Language SBML version 2.
      
reaction ::= 



rule_name ::= 

|   name  


basic_reaction ::= 



reactants ::= 



solution ::= 

| _


enumeration ::= 

| _


kinetics ::= 






  
Example.

,



biocham: a+b+c/d,e => f.






biocham: _ =[a]=> b.






biocham: list_reactions.


[0] MA(1) for (a+b+c)/(d,e)=>f
[1] MA(1) for a=>a+b





biocham: list_ode.


[0] d(b)/dt=a-a*b*c
[1] d(f)/dt=a*b*c
[2] d(c)/dt= - (a*b*c)
[3] d(a)/dt= - (a*b*c)
[4] d(d)/dt=0
[5] d(e)/dt=0
[6] init(b=0)
[7] init(f=0)
[8] init(c=0)
[9] init(a=0)
[10] init(d=0)
[11] init(e=0)












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.2. 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.



3.3. 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 all the reaction rules that match a name prefix.
  



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 kinetics the product of the kinetics (using mass action law kinetics if MA, MAI, MM or Hill kinetics).
  



list_reactions.


lists the current set of reaction rules.



list_rules.


lists the current set of reaction and influence rules.



list_reactions_with_species({object1, ...,
        objectn}).


Lists reactions having a reactant, product or inhibitor in the set of molecular species Objectset.



list_reactions_with_reactant({object1, ...,
        objectn}).


Lists reactions having a reactant in the set of molecular species Objectset.



list_reactions_with_product({object1, ...,
        objectn}).


Lists reactions having a product in the set of molecular species Objectset.



list_reactions_with_catalyst({object1, ...,
        objectn}).


Lists reactions with a catalyst (i.e. both reactant and product) in the set of molecular species Objectset.



list_reactions_with_strict_catalyst({object1, ...,
        objectn}).


Lists reactions with a strict catalyst (i.e. same stoichiometry as reactant and product) in the set of molecular species Objectset.



list_reactions_with_autocatalyst({object1, ...,
        objectn}).


Lists reactions with an autocatalyst (i.e. different stoichiometry as reactant and product) in the set of molecular species Objectset.



list_reactions_with_inhibitor({object1, ...,
        objectn}).


Lists reactions with an inhibitor in the set of molecular species Objectset.



3.4. Reaction graph visualization


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.
  



The graph is either drawn from left to right (default) or from top to bottom.


option(left_to_right:yes).


draw_graph.


Draws the current graph.






Options.




left_to_right: yesno


Draws the graph from left to
    right (default) instead of top to bottom, with species present in the initial state first.







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.
  



3.5. Graph editor


Different graphs can be created from Biocham models and manipulated and more importantly visualized using the third-party Graphviz visualization tool.
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.



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 ::= 



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=a*k/(1+b)
[1] d(a)/dt=0
[2] init(b=0)
[3] init(a=1)
[4] par(k=1)





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 graphs


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



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].
influence_graph.


builds the influence graph between molecular species of the current model without distinguishing between reaction and influence rules.



draw_influences.



    draws the influence graph between molecular species of the current model. Equivalent to influence_graph. draw_graph.
    



  
Example.

  


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






biocham: draw_influences.


Simulation results




    






  



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).


option(force_graph:no).


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







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.



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.4,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(Expression: arithmetic_expression, parameter_name1 = arithmetic_expression1,
       ...,
       parameter_namen =
       arithmetic_expressionn).


Introduce a special kind of event that will be fired when the time
    reaches the given threshold. They may be used for efficiency reason during
    numerical integration.



list_events.


lists all the declared events.



delete_events.


deletes 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/1 ⁄ load_sbml/1
    ⁄ load_model_from_ode/1 ⁄ load_table/1,
    depending on the file extension
    (respectively .bc, .xml, .ode, .csv
    – assuming no extension is .bc).
    If you need to just import an ode system, use import_ode/1.






Options.




use_sbml_id: yesno


Use the sbml_id for the import of all sbml object instead of their names







add(input_file).



    acts as the corresponding add_biocham/1 ⁄ add_sbml/1
    ⁄ add_model_from_ode/1 ⁄ add_table/1,
    depending on the file extension
    (respectively .bc, .xml, .ode, .csv
    – assuming no extension is .bc).






Options.




use_sbml_id: yesno


Use the sbml_id for the import of all sbml object instead of their names







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.



new_model(ModelName).


opens a new fresh model with a selected name.



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.



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.



check_model.


Checks whether the current model is well formed and strict, (see TCS15_FGS).If it is not the case, proposes a witness.



check_polynomiality.


Check whether the current ODE system is polynomial.



correct_model.


Improve wellformedness and strictness of a CRN model by removing reactions of null kinetics, splitting reactions that are implicitely reversible in two and adding correct annotation for modifiers.



6.2. SBML and SBML-qual files




option(use_sbml_id:no).


load_sbml(input_file).



    acts as load_biocham/1, this import the species, reactions, parameters,
    compartments, events and initial states from the provided SBML .xml file. By default,
    we use the sbml Names of the species as biocham identifier. When this leads to naming
    conflicts, we suffix these names with the sbml Id.  The option use_sbml_id (default:
    no) allows to bypass this naming and simply use the id as identifier.
    Rq: Notes and annotations are imported but are not accessible for the user.
  






Options.




use_sbml_id: yesno


Use the sbml_id for the import of all sbml object instead of their names







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 Equation Models


Biocham can also manipulate ODE systems, import and export ODE systems in XPP syntax, and export them in LaTeX.

Remark.  XPP format has restrictions on names and does not distinguish between upper case and lower case letters. For that reason, and in order to avoid misinterpretations, when an XPP file is imported, all names are read in lower case.

The ODE simulation of a Biocham model proceeds by creating an ODE system if there is none, and deleting it after the simulation.
It is worth noting that if there is an ODE system present, it is the current ODE system that is simulated, not the Biocham model.
Biocham can 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.
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 or ODE system.
    



  
Example.

  


biocham: a=>b.






biocham: list_ode.


[0] d(b)/dt=a
[1] d(a)/dt= -a
[2] init(b=0)
[3] init(a=0)




    






  



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*a^2*b*k
[1] d(a)/dt= - (2*a^2*b*k)
[2] d(b)/dt= - (a^2*b*k)
[3] init(c=0)
[4] init(a=0)
[5] init(b=0)
[6] par(k=10)





biocham: export_ode('test2.ode').






biocham: load_reactions_from_ode('test2.ode').






biocham: list_model.


MA(k) for b+2*a=>3*c.
parameter(
  k = 10
).





biocham: load_influences_from_ode('test2.ode').






biocham: list_model.


3*(a^2*b*k) for a,b -> c.
1*(a^2*b*k) for a,b -< b.
2*(a^2*b*k) for a,b -< a.
initial_state(c=0).
initial_state(a=0).
initial_state(b=0).
parameter(
  k = 10
).





biocham: list_ode.


[0] d(a)/dt= - (2*a^2*b*k)
[1] d(b)/dt= - (a^2*b*k)
[2] d(c)/dt=3*a^2*b*k
[3] init(a=0)
[4] init(b=0)
[5] init(c=0)
[6] par(k=10)














option(import_reactions_with_inhibitors:yes).


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 (default yes).







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.



infer_hidden_molecules: yesno


Run infer_hidden_molecules (1-X => X_m) before loading 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




option(second_order_closure:no).


ode_system.



    builds the set of ODEs of the current model.
  






Options.




second_order_closure: yesno


Compute (normal) second order moment closure ODEs. Introduce covariance
      variables and functions to visualize plus⁄minus standard deviation







import_ode(input_file).


imports a set of ODEs.



new_ode_system.


creates an ODE system.



delete_ode_system(name).


deletes an ODE system.



clear_ode_systems.


Remove all existing ODE systems.



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.





option(infer_hidden_molecules:no).


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.



infer_hidden_molecules: yesno


Run infer_hidden_molecules (1-X => X_m) before loading reactions







load_influences_from_ode_system.


Replaces the current reactions with those from
    add_influences_from_ode_system/0.



remove_fraction.


Remove the rational fraction in the current ODE system by multiplying the derivative
  by the GCD. WARNING: the resulting ODE system is NOT equivalent to the starting one.



infer_hidden_molecules.


Tries to infer hidden molecules eliminated by linear invariants.



ode_parameter(name1 = number1,
       ...,
       namen =
       numbern).


sets the value of a parameter in the current set of ODEs.



ode_function(name1 = arithmetic_expression1,
       ...,
       namen =
       arithmetic_expressionn).


defines a (parameterless) function in the current set of ODEs.



load_model_from_ode(input_file).


Clean the current model and import the .ode⁄.xpp file then call add_model_from_ode/1.



add_model_from_ode(input_file).


Import the .ode⁄.xpp file and infer the corresponding reactions or influences system.



export_xpp(output_file).


Exports the ODE system of the current model as an XPP file.



read_xpp(input_file).


Imports an ODE system from an XPP file.



Part II: Qualitative Analysis and Synthesis
Chapter 7
Static Analyses
7.1. Graph properties



      Analysis of some graph properties of the current model.
      
pathway(object1, object2).


Gives one reaction pathway from object1 to object2 if one exists in the directed reaction graph of the current model (for more complex queries, see next section on Computation Tree Logic model-checking).



list_input_species.


Lists the molecular species that are neither a reaction product apart from strict catalysts, nor the target of an influence rule, in the curent model. Reactants that are not catalysts, and strict catalysts, i.e. catalysts with same stoichiometry in left-hand and right side of reactions, are allowed as inputs.



list_strict_input_species.


Lists the molecular species that are strict catalysts, i.e. neither a strict product, nor strict reactant, nor the target of an influence rule in the curent model.



list_source_species.


Lists the molecular species that are neither a reaction product nor the target of an influence rule in the curent model.



list_sink_species.


Lists the molecular species that are neither a reactant nor a positive source of an influence rule in the curent model.



7.2. Dimensional analysis


unit_type ::= 

| time
| volume
| substance


set_units(Type: unit_type, Unit: name).


Set the default units for the given type (time, volume or substance) in the current model.



list_units.


Lists the default units of the current model.




   Dimensional analysis infers dimensions for model parameters and checks their consistency.
   In Biocham, only time and volume dimensions are considered.
   The dimension of a molecular concentration is volume$^{-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 according to ODEs.
     An error is raised if some parameters appear in expressions with inconsistent dimensions.
     This can happen in a model if some  multiplicative parameter with value equal to 1 is omitted for example.
    



  
Example.

    


biocham: MM(v, k) for A => B.






biocham: parameter(k = 2, v = 3).






biocham: list_dimensions.


v has dimension time^(-1).volume^(-1)
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.3. 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.







search_efms.


computes and displays the T-invariants of the network up to the maximal
    flux given by the option conservation_size and defaulting to 4. Such
    T-invariants can be seen as a way to compute Extreme Rays of the cone of
    Elementary Flux Modes.






Options.




conservation_size: integer


Set the maximal stoichiometric size for a P-invariant.







7.4. 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).
| not_species
Search reduction while forbidding to merge species
together.
| 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).




option(show_support: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)



show_support: yesno


determine if the printing describe the SEPI or only its
  support.



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 A+B=>A+C.





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
{E+S => ES} -> {A+B => A+C}
{ES => E+S} -> deleted
{ES => E+P} -> {A+B => A+C}
Number of reductions: 1









7.5. 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 ligand 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.6. 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~{p1}' = Cdc2~{p1}
Cyclin' = ε^(- 4) * Cyclin
Cdc2-Cyclin~{p1,p2}' = Cdc2-Cyclin~{p1,p2}
Cdc2-Cyclin~{p1}' = ε^(- 2) * Cdc2-Cyclin~{p1}
Cyclin~{p1}' = ε^(- 2) * Cyclin~{p1}


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


found a complete equilibration leading to the rescaling:
Cdc2' = ε^(- 4) * Cdc2
Cdc2~{p1}' = ε^(- 2) * Cdc2~{p1}
Cyclin' = ε^(- 2) * Cyclin
Cdc2-Cyclin~{p1,p2}' = Cdc2-Cyclin~{p1,p2}
Cdc2-Cyclin~{p1}' = ε^(- 2) * Cdc2-Cyclin~{p1}
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.7. 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 may be non-degenerate multistationarity, positive circuit detected.









check_oscillations.


Checks a necessary condition for the periodic behaviour of the current model.



7.8. Rate independence


Under the condition that a CRN is well-formed (i.e. the of species occurring in the rate of a reaction is the set of its reactants, catalysts and inhibitors [FGS15tcs], Biocham can check sufficient conditions that ensure that the outputs of a CRN are independent of the reaction rates [DFS20cmsb].
test_rate_independence.


Tests graphical sufficient conditions for rate independence of the current model for all output species (assuming well-formed kinetics).



test_rate_independence_invariants.


Tests invariant sufficient conditions for rate independence of the current model for all output species (assuming well-formed kinetics).



7.9. Mean stochastic behavior


Biocham can check sufficient conditions to ensure that, for some molecular species of the current reaction model, the mean stochastic trace is given by the ODE simulation trace, at all time points, and for any conversion factor (i.e. with small molecule numbers, unlike Kurtz's limit theorem). The command returns the list of molecular species for which that property is guaranteed to hold. The condition basically checks that those species are produced by reactions with mass action law kinetics and that polymolecular reactions are restricted to catalytic synthesis reactions with disjoint sets of ancestor species.
test_stochastic_mean_ode.


Lists the molecular species of the current reaction model for which the ODE simulation trace is guaranteed to be equal to the mean stochastic trace at all time points, for any conversion factor. This is always the case of the strict input species (see list_strict_input_species/0) which are listed apart.



Chapter 8
Boolean Dynamics, Verification and Synthesis
8.1. Attractors in the ground Boolean state transition graph


The following commands refer to the Boolean dynamics of a BIOCHAM model, possibly combining reaction rules and influence rules.
      The Boolean semantics 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].
The notion of Boolean states considered in this first section is the one of complete ground states defined by the presence or absence of each molecular species,
unlike the notion of symbolic representation of partial states represented by Boolean constraints and considered in the following sections.
list_stable_states.



    lists stable steady ground Boolean states of the ground state transition graph corresponding
    to the Boolean semantics of the current model.
  






Options.




boolean_semantics: boolean_semantics


Use positive or negative boolean semantics for inhibitors.



boolean_state_display: boolean_state_display


choice of display of the boolean states.







list_tscc_candidates.



    lists possible representative states of Terminal Strongly Connected
    Components (TSCC) of the ground state transition graph corresponding to the
    positive semantics of the current model.
  






Options.




boolean_state_display: boolean_state_display


choice of display of the boolean states.







8.2. Temporal logic properties of the symbolic Boolean state transition graph


The Computation Tree Logic CTL can be used to express qualitative properties of the Boolean dynamics of BIOCHAM model in a given (set of) initial states [CCDFS04tcs].
CTL extends propositional logic, used to represent set of states symbolically by Boolean constraints, with modal operators,
       used 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)) expresses the reachability of a stable set of states where protein p will always remain present.
It is worth noting that in this setting, a propositional formula provides a symbolic representation for a set of ground Boolean states,
      also called a partial state, or just a state by abuse of notation when there is no ambiguity with the ground Boolean states considered in the previous Section.
      CTL is well suited to analyze attractors, their reachability and transient properties  [TFFT16bi].
It is worth noting that reachability properties in temporal logics are relative to a given set of initial states, they are purely factual and not necessarily causal.
The syntax of CTL formulas, plus 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  )
| must_reach(  ctl  )
| steady(  ctl  )
| stable(  ctl  )
| checkpoint(  ctl  ,  ctl  )
| checkpoint2(  ctl  ,  ctl  )
| oscil2(  ctl  )
| oscil3(  ctl  )
| oscil(  ctl  )
| current_spec


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 not(g)⁄\ 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.



clear_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).






Options.




boolean_initial_states: boolean_initial_states


specifies whether the truth value of a formula is for all or some completion of the initial states (present⁄absent⁄undefined).







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)









8.3. Verification of CTL and LTL formulae


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).


The negative Boolean semantics interpret reaction inhibitors with a negation in the condition. The positive Boolean semantics ignore reaction inhibitors.
boolean_initial_states ::= 

| all
| some
| all_then_some
| present



  The default value is all. The value
  all_then_some tries for all initial states, and if that fails, for
  some. The value present treats undetermined species as absent.
boolean_state_display ::= 

| present
| table
| vector


Display boolean states as a present command. (default), or a TAB separated table of TRUE⁄FALSE values, or an association-list with 0⁄1 values.


option(boolean_initial_states:all).


CTL formulae can be checked in all or some initial states.


option(boolean_trace:no).


Display either a Boolean trace proving the formula in some initial state, or a counter example initial state falsifying the formula (default value is no trace).


option(nusmv_topological_order:no).


Ordering of variables for the NuSMV model-checker.


option(query:current_spec).




option(boolean_state_display:present).


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 boolean_initial_states option.
  






Options.




query: ctl


Query to evaluate instead of the current CTL
    specification.



boolean_initial_states: boolean_initial_states


specifies whether the truth value of a formula is for all or some completion of the initial states (present⁄absent⁄undefined).



boolean_trace: yesno


shows a proof trace or a counter example initial state if available  (default value is
      "no").



boolean_semantics: boolean_semantics


choice of positive⁄negative boolean semantics (positive: reaction inhibitors are ignored; negative: the presence of one reaction inhibitor prevents the reaction to proceed.



nusmv_topological_order: yesno


tells NuSMV to keep the order of the variables of the BIOCHAM model for building its internal Binary Decision Diagram.



boolean_state_display: boolean_state_display


choice of display of the boolean states.










  
Example.




biocham: present(a).






biocham: absent(b).






biocham: a=>b.






biocham: a+b=>a.






biocham: check_ctl(query: EX(not(a) \/ EG(not(b))), boolean_trace:yes).


EX(not a\/EG(not b)) is true





biocham: check_ctl(query: EG(not(b)), boolean_trace:yes).


Trace:
present({a}).

EG(not b) is false





biocham: check_ctl(query:reachable(b),boolean_trace: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









A Linear Time Logic (LTL) formula contains no path quantifier. The syntax of LTL formulae is given in the next chapter for the more general setting of First-Order LTL (FO-LTL) formulae with numerical constraints.
      Here a boolean LTL formula is true if it is satisfied by all paths and all boolean initial states.
      This can be changed by double negation with the option boolean_initial_states: some to verify that there exists an initial state and a path satisfying the formula.
      Unlike CTL, the path quantifier cannot be decoupled from the quantifier on the initial state.
check_ltl(Query: foltl).



    evaluates the given boolean LTL specification on the current model by
    calling the NuSMV model-checker in Bounded Model-Checking mode.
  






Options.




boolean_initial_states: boolean_initial_states


specifies whether the truth value of a LTL formula is for all or some completion of the initial states (present⁄absent⁄undefined) and for all or some paths.



boolean_trace: yesno


shows a proof trace or a counter example initial state if available  (default value is
      "no").



boolean_semantics: boolean_semantics


Use positive or negative boolean semantics for inhibitors.



boolean_state_display: boolean_state_display


choice of display of the boolean states.










  
Example.




biocham: present(a).






biocham: absent(b).






biocham: a=>b.






biocham: a+b=>c.






biocham: check_ltl(reachable(c), boolean_trace:yes).


Trace:
present({a}).
present({b}).
present({b}).

Reactions used for that path:
MA(1)for a=>b
reachable(c) is false





biocham: check_ltl(reachable(c), boolean_initial_states: some, boolean_trace:yes).


Trace:
present({a, c}).

reachable(c) 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.4. 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 RAFPH+'RAF~{p1}'=>'RAF~{p1}-RAFPH'
MA(1) for 'RAF~{p1}-RAFPH'=>RAFPH+'RAF~{p1}'
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'=>MEKPH+'MEK~{p1}'
MA(1) for 'MAPK~{p1}-MAPKPH'=>MAPK+MAPKPH
MA(1) for 'MAPK~{p1,p2}-MAPKPH'=>MAPKPH+'MAPK~{p1}'





biocham: dose_response('RAFK',1e-5,1e-3, time:200, show:'MAPK~{p1,p2}').


Simulation results












8.5. 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.6. 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).


% Maximum K used: 1
% minimum number of samples for h=1: 18

% 102 samples (max h ~ 5.666666666666667)
P -< P

% 0 samples (max h ~ 0)
% _ -> P

% 30 samples (max h ~ 1.6666666666666667)
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.
initial_state(R=1).
initial_state(P=1).
parameter(
  k1 = 2,
  k2 = 2,
  k3 = 1
).





biocham: pac_learning(library:examples/circadian_cycle/bernot_comet.bc, 100, 10).


% Maximum K used: 2
% minimum number of samples for h=1: 54

% 196 samples (max h ~ 3.6296296296296298)
G,PC -< G

% 29 samples (max h ~ 0.5370370370370371)
/ PC,G -> G

% 456 samples (max h ~ 8.444444444444445)
L -< L

% 134 samples (max h ~ 2.4814814814814814)
/ L -> L

% 176 samples (max h ~ 3.259259259259259)
L,PC -< PC
PC / G -< PC

% 9 samples (max h ~ 0.16666666666666666)
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


polynomial ::= 



monomial ::= 

|   name  



  Biocham is interfaced to 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 is the default one.

 The ODE simulation of a Biocham model proceeds by creating an ODE system if there is none, and deleting it after the simulation.
It is worth noting that if there is an ODE system already present (e.g. created by import_ode), it is the current ODE system that is simulated, not the Biocham model.
The stochastic simulation of a Biocham model is specific to reaction and influence models which can be interpreted by a Continuous Time Markov Chain (CTMC).
      The stochastic simulation algorithm cannot be directly used on an ODE system, but on an equivalent reaction system (e.g. automatically inferred with command load_reactions_from_ode/1).
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
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
  the initial concentrations are converted in molecule numbers by using the stochastic\_conversion option (default 100) and 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
Experimental implementation of Rosenbrock method (not compatible with search_parameters and robustness and sensitivity commands).


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(steps:1).




option(stochastic_conversion:100).




option(stochastic_bound:1000000.0).




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



steps: integer


number of steps at which the output is evenly sampled



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 (default 100).



stochastic_bound: number


Maximum number of molecules of one molecular species allowed in stochastic simulations (default 1e6).



stochastic_thresholding: number


Do not write (but still compute) stochastic steps below one fraction of
    the total time (default 1000)



filter: filter


filtering function for the trace







continue.


Continues the numerical simulation by extending the time horizon by option "time" units.






Options.




time: time


time to add to the numerical integration







9.2. Hybrid simulations


This section describes a hybrid numerical simulation
  method, defined either statically by distinguishing the continuous reactions from the stochastic reactions,
       or dynamically using a strategy based on the number of particles and velocity of the reactions.
            Species appearing in both continous and stochastic reactions are
  called hybrid species. The total particle count of one species is 
  calculated and plotted in the simulation trace.


  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 initial_state(Variable = Value)) 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 initial_state(_,_)) 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(plot_table:'').




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.




plot_table: name


Determined the table to be plotted



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_svg(output_file).


plots the current trace in a SVG file



export_plot_to_pdf(output_file).


plots the current trace in a PDF 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 by slowly varying the dose,
under the hypothesis that the input dose is a catalyst not affected itself 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.1-1.0e-7)/6 for _=>a.
variation_b--b*log((1.0e-20+0.1)/(1.0e-20+1.0e-7))/6 for b=>2*b.
initial_state(a=1.0e-7).
initial_state(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 (the dose) of the input object which must be a catalyst not affected by the system, 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 for all values of the dose. 
     The slow variation of the dose is performed over a time horizon of 100 fold that simulation time.






Options.




time: time


time horizon necessary to approximate the steady state response; the variation of the dose is performed over a time horizon of 100 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 which must be a catalyst not affected by the system, 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 slow increasing and decreasing variations of the dose are performed over a time horizon of 100 fold the simulation time.
     Unlike a real bifurcation diagram that draws the stable zeros of the differential function of the response, the steady states are not computed here but approximated by simulation.






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.



set_initial_last_state.


sets the last state as new initial state.



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 properties or 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].
FO-LTL(Rlin) formulas are evaluated on infinite traces obtained from finite traces by automatically completing them with a loop on the last state. Finite numerical data time series can be generated by simulation or loaded from biological experiment data.
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). Verifying the satisfiability of an FO-LTL(Rlin) formula thus amounts to compute the validity domain of its variables.
The symbol == tests equality of real numbers on a discrete trace by interpolation and should be prefered to exact equality = as explained and illustrated in the example below.
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 ::= 

|   name  


time_value ::= 



time_value_list ::= 

| 
      [ time_value
       , ... ,
       time_value ] 


parameter_between_interval ::= 



variable_objective ::= 

|   name  ->  number  


Equality over the real numbers is the equality of floating point numbers. A predicate for equality by interpolation with next time point, noted A==B, is also defined on a discrete trace by A = B \⁄ (A > B ⁄\ X(A < B)) \⁄ (A < B ⁄\ X(A > B))
      i.e. the interpolation equality holds in any state satisfying just one inequality with the opposite inequality in its next state (Rolle Theorem). The interpolation equality should be used for instance when fixing the Time state variable to a value that is not present in the trace (see example below).
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.
The foltl_precision option (default 12, i.e. $10^{-12}$) is used for converting and solving FO-LTL constraints over the integers instzead of floating point numbers.


option(foltl_precision: 12).




option(foltl_magnitude: 5).


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: MA(k) for A => B.






biocham: present(A,a).






biocham: parameter(a=1, k=1).






biocham: option(time:5).






biocham: numerical_simulation.






biocham: plot.


Simulation results





biocham: validity_domain(F(Time==1 /\ A=x /\ B=y)).


x=0.371929/\y=0.628071





biocham: validity_domain(F(Time=1 /\ A=x /\ B=y)).


false





biocham: validity_domain(F(Time=T /\ A==B)).


T=0.689051





biocham: validity_domain(F(Time=T /\ A=B)).


false





biocham: validity_domain(F(A==B /\ A=x /\ B=y)).


x=0.502052/\y=0.497948











  






Options.




foltl_precision: number


precision for conversion to integers



foltl_magnitude: number


order of magnitude for << and >> operators







satisfaction_degree(Formula: foltl, [variable_objective1, ...,
        variable_objectiven]).



    computes the continuous satisfaction degree in the interval [0,+∞) of an
    FOLTL(Rlin) formula on the current trace with respect to objective values given for each free variable of the formula (variables without objective should be existentially quantified in the formula). The satisfaction 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).
     



  
Example.




biocham: satisfaction_degree(F(Time==T /\ A=x /\ B=y), [T -> 0.689, x -> 0.5, y-> 0.5]).


0.997106





biocham: satisfaction_degree(F(Time==T /\ A=x /\ B=y), [T -> 1, x -> 0.5, y-> 0.5]).


0.877357




  






  






Options.




foltl_precision: number


precision for conversion to integers



foltl_magnitude: number


order of magnitude for << and >> operators







ltl_pattern(function_prototype1 = foltl1,
       ...,
       function_prototypen =
       foltln).


defines new F0-LTL formula patterns.



list_ltl_patterns.


lists all known FO-LTL formula patterns.



expand_ltl(Formula: foltl).


shows the expansion in FO-LTL of a formula pattern.



delete_ltl_pattern(functor1,
       ...,
       functorn).


deletes some LTL patterns. Either the arity is given with the pattern name⁄arity, or all LTL patterns with
  the given functor are deleted.



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)))).


10.3. Parameter sensitivity, robustness and parameter optimization w.r.t. FO-LTL(Rlin) properties




option(robustness_coeff_var: 0.1).




option(robustness_samples: 30).




option(positive_sampling: yes).




option(openmpi_procs: 1).


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 in multiple ways: either to search parameter values in high parameter dimension to satisfy a temporal property [RBFS11tcs], or compute robustness measures and parameter sensitivity indices with respect to parameter perturbations estimated by sampling [FS18cmsb], or to visualize the landscape in two dimensions.
Those commands involve numerical simulations which are executed by samplings parameter values according to a normal distribution law around their current values with a given coefficient of variation (0.1 by default), and taking absolute values (option positive_sampling by default). The number of parameter set samples to consider is given by the option robustness_samples (default 30) which will be multiplied by the number of parameters. You can use the seed/1 command to try different samplings and check whether the results are similar, otherwise the number of samples should be increased.
The simulations are automatically parallelized by OpenMPI. The number of processors used is given by an option (by default 0 for automatic, temporarily replaced by 1 by default because of a bug on Ubuntu)
robustness(Formula: foltl, [parameter_name1, ...,
        parameter_namen], [variable_objective1, ...,
        variable_objectiven]).



    computes the robustness degree defined as the mean satisfaction degree (truncated to at most 1 for each parameter set sample simulation) [RBFS11tcs],
      of an FOLTL formula Formula, 
     given with a list of [variable_objective1, ...,
        variable_objectiven] values for the free variables of the formula (variables without any objective value should be existentially quantified)
    with respect to a variation of the 
     [parameter_name1, ...,
        parameter_namen] (i.e. extrinsic variability).The list of parameters can also be empty to measure the robustness with respect to stochastic simulations (i.e. intrinsic variability).
    The robustness is  estimated by sampling with a normal perturbation law 
       around the current parameter values according to the coefficient of variation option (default 0.1 i.e. 10 per cent).




  
Example.

     Continuing the previous example:



biocham: robustness(F(A==B /\ A=x /\ B=y), [a, k], [x -> 0.5, y-> 0.5]).


 On [a,k] robustness degree: 0.951512 sensitivity: 0.048488  (computation time 4.286000 s)











  






Options.




time: time


time horizon of the numerical integration



method: method


method for the numerical solver



robustness_coeff_var: number


coefficient of variation (default 0.1, i.e. 10%)



positive_sampling: yesno


taking the absolute values of the sampled values (default yes)



robustness_samples: integer


number of samples (default 30) per parameter



openmpi_procs: integer


Number of processors to use with OpenMPI.



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







sensitivity(Formula: foltl, [parameter_name1, ...,
        parameter_namen], [variable_objective1, ...,
        variable_objectiven]).



    computes sensitivity indices for each parameter in the list [parameter_name1, ...,
        parameter_namen],
     according to a normal distribution around its current value with some given coefficient of variation (default 0.1),
      as the mean satisfaction degree of the formula Formula
     with [variable_objective1, ...,
        variable_objectiven] objective values for its free variables.
     This command is just a shorthand for calling robustness repeatedly on each single parameter in the list one at a time (i.e. OAT method based on the mean).






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 parameters (default 0.1)



positive_sampling: yesno


taking the absolute values of the sampled values (default yes)



robustness_samples: integer


number of samples (default 30) per parameter



openmpi_procs: integer


Number of processors to use with OpenMPI.



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







sensitivity(Formula: foltl, [variable_objective1, ...,
        variable_objectiven]).



    computes the sensitivity indices of all model parameters. This command is a shorthand for calling the previous command with all parameters.
     



  
Example.

     The command computes sensitivity indices for both parameter a and k with default coefficient of variation of 0.1:



biocham: sensitivity(F(A==B /\ A=x /\ B=y), [x -> 0.5, y-> 0.5]).


 On [a] robustness degree: 0.953328 sensitivity: 0.046672  (computation time 3.848000 s)
 On [k] robustness degree: 0.985727 sensitivity: 0.014273  (computation time 3.840000 s)
[a,k] by order of decreasing sensitivity.


















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 (default 0.1)



positive_sampling: yesno


taking the absolute values of the sampled values (default yes)



robustness_samples: integer


number of samples (default 30) per parameter



openmpi_procs: integer


Number of processors to use with OpenMPI.



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.

     This command allows us to tune the kinetic parameter k to obtain a desired timing in the previous example:



biocham: search_parameters(F(Time==T /\ A==B), [0 <= k <= 2], [T -> 1]).


Time: 4.065 s
Stopping reason: Fitness: function value -1.08e-02 <= stopFitness (1.00e-04)
Best satisfaction degree: 1.010949
[0] parameter(k=0.6833050831000044)











  






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)



openmpi_procs: integer


Number of processors to use with OpenMPI.







search_parameters([parameter_between_interval1, ...,
        parameter_between_intervaln], [search_condition1, ...,
        search_conditionn]).


similar to search_parameters/3 but using a list of multiple 
    [search_condition1, ...,
        search_conditionn], each condition is given as a triple composed of an FOLTL formula, a list of particular parameters instantiations (e.g. describing a mutant), and  objective values for the
    variables of the FOLTL 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.058 s
Stopping reason: Fitness: function value 4.68e-06 <= stopFitness (1.00e-04)
Best satisfaction degree: 0.999995
[0] parameter(ka=1.5000002341053575)















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)



openmpi_procs: integer


Number of processors to use with OpenMPI.









option(resolution: 10).


scan_parameters(Formula: foltl, Parameter1: parameter_between_interval, Parameter2: parameter_between_interval, [variable_objective1, ...,
        variable_objectiven]).



    returns some statistics and draws the landscape of the satisfaction degrees (truncated to 1) of a FO-LTL(Rlin) formula obtained by varying two parameters in intervals given, possibly on loagarithmic scale. Note that the number of calls to numerical integration is the square of the resolution option value N. By default the image and the data are saved in files landscape_Parameter1_Parameter2_N.png and .txt respectively. Note that the drawing shows an overapproximation of the best values found for the parameters. The precise values are available in file  landscape_Parameter1_Parameter2_N.txt.
 



  
Example.

 


biocham: MA(k)for a => b. present(a,10).






biocham: MA(l) for b => a.






biocham: scan_parameters(F(Time > T /\ a > b), (0 <= k <= 2), (0 <= l <= 2), [T -> 5], resolution:3).


Time: 28.999 s
Mean satisfaction (robustness): 0.380110, standard deviation: 0.584952
Maximum satisfaction: 1.000000
Best solutions found with 0.000000<=k<=1.000000, 0.000000<=l<=2.000000
Simulation results




 






     






Options.




logscale: axes


parameter sampling and axes on a logarithmic scale



resolution: number


number of values tried for each parameter (default is 10)



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 and minimum value on logarithmic scale



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



openmpi_procs: integer


Number of processors to use with OpenMPI.







scan_parameters(Formula: foltl, Parameter1: parameter_between_interval, Parameter2: parameter_between_interval, [variable_objective1, ...,
        variable_objectiven], Filename: output_file).



    specifies in addition to save 
     the image in file filename.png and the data in file filename.txt.
     






Options.




logscale: axes


sampling and axes on a logarithmic scale



resolution: number


number of values tried for each parameter (default is 10)



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



openmpi_procs: integer


Number of processors to use with OpenMPI.







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_n*z_p
[1] d(z_n)/dt=x_n+y_n-z_n-fast*z_n*z_p
[2] d(x_p)/dt=0
[3] d(y_p)/dt=0
[4] d(x_n)/dt=0
[5] d(y_n)/dt=0
[6] d(x)/dt=0
[7] d(y)/dt=0
[8] d(temp0_p)/dt=0
[9] d(temp0_n)/dt=0
[10] d(temp1_p)/dt=0
[11] d(temp1_n)/dt=0
[12] init(z_p=0)
[13] init(z_n=0)
[14] init(x_p=1)
[15] init(y_p=3)
[16] init(x_n=0)
[17] init(y_n=0)
[18] init(x=1)
[19] init(y=3)
[20] init(temp0_p=0)
[21] init(temp0_n=0)
[22] init(temp1_p=0)
[23] init(temp1_n=0)
[24] par(fast=1000)





biocham: list_model.


MA(1) for x_p=>x_p+z_p.
MA(1) for y_p=>y_p+z_p.
MA(1) for z_p=>_.
MA(1) for x_n=>x_n+z_n.
MA(1) for y_n=>y_n+z_n.
MA(1) for z_n=>_.
MA(fast) for z_n+z_p=>_.
initial_state(x=1).
initial_state(y=3).
initial_state(temp0_p=0).
initial_state(temp0_n=0).
initial_state(temp1_p=0).
initial_state(temp1_n=0).
initial_state(x_p=1).
initial_state(x_n=0).
initial_state(y_p=3).
initial_state(y_n=0).
initial_state(z_p=0).
initial_state(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 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 with full generality.
The option for quadratic reduction restricts the synthesis to reactions with at
 most two reactants (the default is not, since this process may be costly).
 Two methods are available to perform this a posteriori reduction: natively
 or using an external SAT solver.


option(quadratic_reduction:fastnSAT).




option(quadratic_bound:carothers).


Another option is the lazy introduction of molecular species for negative real values (the default is yes).


option(lazy_negatives:yes).




option(negation_first: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.




quadratic_reduction: reduction_methods


Determine if the quadratic reduction for synthesizing reactions with at most two reactants has to be performed



quadratic_bound: reduction_bounds


Use the bound from Carother's article or that from Bychkov and Pogudin



lazy_negatives: yesno


Switch between systematic and lazy introduction of molecular species for negative values



negation_first: yesno


Determine if the negation of variables is perform before or after the quadratization







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.




quadratic_reduction: reduction_methods


Determine if the quadratic reduction for synthesizing reactions with at most two reactants has to be performed



quadratic_bound: reduction_bounds


Use the bound from Carother's article or that from Bychkov and Pogudin



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.0) for D=>D+H+output.
MA(1.0) for F=>F+H+output.
MA(1.0) for _=>D+F.
initial_state(output=4).
initial_state(I=4).
parameter(
  fast = 1000
).





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.


MA(fast) for costime_m+costime_p=>_.
MA(fast) for J_m+J_p=>_.
MA(1.0) for J_p=>J_p+costime_p.
MA(1.0) for J_m=>J_m+costime_m.
MA(1.0) for costime_m=>J_p+costime_m.
MA(1.0) for costime_p=>J_m+costime_p.
initial_state(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.


MA(fast) for costime1_m+costime1_p=>_.
MA(fast) for N_m+N_p=>_.
MA(fast) for K_m+K_p=>_.
MA(1.0) for K_p=>K_p+N_p+costime1_p.
MA(1.0) for K_m=>K_m+N_m+costime1_m.
MA(1.0) for N_m=>K_p+N_m.
MA(1.0) for N_p=>K_m+N_p.
initial_state(costime1_p=2).
initial_state(N_p=1).
initial_state(J=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.


MA(fast) for cosx_m+cosx_p=>_.
MA(fast) for K_m+K_p=>_.
MA(1.0) for K_p+x=>K_p+cosx_p+x.
MA(1.0) for K_m+x=>K_m+cosx_m+x.
MA(1.0) for cosx_m+x=>K_p+cosx_m+x.
MA(1.0) for cosx_p+x=>K_m+cosx_p+x.
MA(1.0) for x=>_.
initial_state(cosx_p=1).
initial_state(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 with initial values defined by polynomials of the input. PIVPs are written with the following grammar:
pivp ::= 

|   pivp  ;  pivp  


pode ::= 

| d(  name  )/dt=  polynomial  


compile_from_pivp(PIVP: pivp, name1,
       ...,
       namen).



    creates a CRN that immplements a function of time, specified by the variable "Output" of the solution of the PIVP.
  






Options.




quadratic_reduction: reduction_methods


Determines if the quadratic reduction has to be performed



quadratic_bound: reduction_bounds


Use the bound from Carother's article or that from Bychkov and Pogudin








  



  
Example.
 Compilation of the Hill function of order 2 as a function of time, $h(t)=t^2⁄(1+t^2)$ 




biocham: option(quadratic_reduction:native).






biocham: compile_from_pivp((0.0,d(h)/dt=2*n^2*t;1.0,d(n)/dt= -2*n^2*t;0.0,d(t)/dt=1),[h]).






biocham: list_model.


MA(2.0) for n+nt=>h+nt.
MA(1.0) for n=>n+nt.
MA(2.0) for 2*nt=>nt.
initial_state(n=1.0).
parameter(
  fast = 1000,
  input = 4
).





biocham: numerical_simulation(time:10).






biocham: plot(show:{h}).


Simulation results





biocham: plot(show:{h},logscale:x).


Simulation results




  







compile_from_pivp(PIVP: pivp, Input: name, name1,
       ...,
       namen).



  creates a CRN that implements a function of the variable "Input". The function is specified by the value of the variable "Output" of the solution of a PIVP at time "Input".
     The "Input" variable must not appear in the PIVP. It is created and  initialized with a parameter named "input".
  






Options.




quadratic_reduction: reduction_methods


Determines if the quadratic reduction has to be performed



quadratic_bound: reduction_bounds


Use the bound from Carother's article or that from Bychkov and Pogudin







compile_from_ode(Input: name, name1,
       ...,
       namen).


Compile the current ODE system in a CRN, Input is usually "time" but could be the name of a new species, Output have to be the name of an existing variable.






Options.




quadratic_reduction: reduction_methods


Determines if the quadratic reduction has to be performed



quadratic_bound: reduction_bounds


Use the bound from Carother's article or that from Bychkov and Pogudin



lazy_negatives: yesno


Switch between systematic and lazy introduction of molecular species for negative values








  



  
Example.
 Compilation of the Hill function of order 2 as a function of an input $h(x)=x^2⁄(1+x^2)$.
  This time sat_species is used (minimizing the number of variables as is the case for native) and then
  sat_reactions that minimizes the number of monomes.



biocham: option(quadratic_reduction:sat_species).






biocham: compile_from_pivp((0.0,d(h)/dt=2*n^2*t;1.0,d(n)/dt= -2*n^2*t;0.0,d(t)/dt=1),x,[h]).






biocham: parameter(input=2).






biocham: list_model.


MA(1.0) for x=>_.
MA(2.0) for n+ntx=>h+ntx.
MA(1.0) for nx=>_.
MA(2.0) for ntx+nx=>ntx.
MA(1.0) for nx+x=>ntx+nx+x.
MA(1.0) for ntx=>_.
MA(2.0) for 2*ntx=>ntx.
initial_state(x=input).
initial_state(n=1.0).
initial_state(nx=input).
parameter(
  fast = 1000,
  input = 2
).





biocham: numerical_simulation(time:10).






biocham: plot(show:{h}).


Simulation results





biocham: plot(show:{h},logscale:x).


Simulation results





biocham: compile_from_pivp((0.0,d(h)/dt=2*n^2*t;1.0,d(n)/dt= -2*n^2*t;0.0,d(t)/dt=1),x,[h],quadratic_reduction:sat_reactions).






biocham: list_model.


MA(1.0) for x=>_.
MA(2.0) for n+ntx=>h+ntx.
MA(1.0) for nx=>_.
MA(2.0) for ntx+nx=>ntx.
MA(1.0) for nx+x=>ntx+nx+x.
MA(1.0) for ntx=>_.
MA(2.0) for 2*ntx=>ntx.
initial_state(x=input).
initial_state(n=1.0).
initial_state(nx=input).
parameter(
  fast = 1000,
  input = 1.0
).




  







reduction_methods ::= 

| no
| native
| fast
| fastnSAT
| sat_species
| sat_reactions


reduction_bounds ::= 

| carothers
| bychkov


quadratize_ode(object1,
       ...,
       objectn).


Modify the ODE system in order to make it quadratic while keeping the prescribed variables



polynomial_ODE.


Check if the current ODE system is polynomial



compile_function(Formula, Input: object).


Compile the desired formula into a CRN by determining a corrsponding IVP, then PIVP and finally a CRN. Typical usage looks like:
  - compile_function(Y=f(X), X). Which set up a CRN such that when parameter(input=X), the final value of Y is f(X). and
  - compile_function(Y=f(Time), Time). Which set up a function of time.






Options.




quadratic_reduction: reduction_methods


Determine if the quadratic reduction for synthesizing reactions with at most two reactants has to be performed







polynomialize_ode.


Starting from the current ODE system, builds an equivalent polynomial ODE system by introducing more variables that compute the non-polynomial terms.




  



  
Example.
 Polynomialization of an ODE system:
 



biocham: d('A')/dt=exp('A').






biocham: polynomialize_ode.


B = exp(A)





biocham: list_ode.


[0] d(B)/dt=B^2
[1] d(A)/dt=B
[2] init(B=1.0)
[3] init(A=0)




  







stabilize_expression(Expression: arithmetic_expression, Output: name, Point).



    Construct a CRN such that the Output species will adjust its concentration
    to satisfy the algebraic Expression in the vicinity of Point.
    




  



  
Example.
 Compilation of a CRN to stabilize the Bring radical:



biocham: stabilize_expression(y^5+y-x,y,[x=2,y=1]).






biocham: list_model.


MA(1.0) for x=>x+y.
MA(1.0) for y=>_.
MA(1.0) for y+y4=>y4.
MA(1.0) for 2*y=>y2+2*y.
MA(1.0) for y2=>_.
MA(1.0) for 2*y2=>y4+2*y2.
MA(1.0) for y4=>_.
initial_state(x=2).
initial_state(y=1).
initial_state(y2=1).
initial_state(y4=1).
parameter(
  fast = 1000,
  input = 1.0
).





biocham: option(time:100,show:{y}).






biocham: dose_response(x,0,4).


Simulation results




  






  
11.3. 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 f+x2=>f+x1+x2.
fast*[x1] for x1=>_.
MA(1) for x1=>x0+x1.
MA(1) for x0=>f+x0.
initial_state(x2= -1).
initial_state(f=1).
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: transfer_function, U: object, Y: object).


compile a transfer function in variable s into a chemical filter.



compile_transfer_function([number1, ...,
        numbern], [number1, ...,
        numbern], 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




      






      

Part IV: Index of Commands
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