Modeling space in BIOCHAM

April - May 2005
Sakina-Dorothée Ayata,
Internship in Contraintes Project, INRIA Rocquencourt, France

During my internship in the Contraintes team (april-may 2005), I worked on space modeling in the Biochemical Abstract Machine BIOCHAM, from the 2.2 and 2.3 versions (the 2.3 version will soon be available). What I have done won't be included in the next release, but will be a base for a further development of space modeling in BIOCHAM.
Those pages aim at presenting my work in a simple way, so you wont need to know how in details, how BIOCHAM works to understand it.

Why modeling space in BIOCHAM ?
Brief presentation of BIOCHAM and context, interest of space modeling.
Principal modifications of BIOCHAM
Presentation of the changes I have introduced in:
- new syntax and molecule localisation
- volume ratio and transport rules recognition
Examples of space modeling with BIOCHAM
Two examples to illustrate space modeling in BIOCHAM:
- p53/Mdm2 network: DNA repair in stress conditions
- Delta/Notch signaling: cellular differentiation
A possible development for the future
Further developments that could be added

Introduction: why modelising space in BIOCHAM ?

Currently with BIOCHAM it is possible to model biochemical reactions, to simulate cellular pathways (with boolean or numerical models) and to ask for temporal properties. Numerical simulations of biochemical networks allow to predict concentration evolutions over time.
Until now, BIOCHAM didn't consider the spatial distribution of the molecules. During simulations, all molecules are mixed together and reactions occur in a molecular soup. But in the cells, molecules are localised and don't play the same roles in different locations. Moreover transport phenomena can induce reaction delays and this must be taken into account for temporal simulations.
The aim of my internship was to introduce space modeling into BIOCHAM and to supply a few examples of its importance. You can find the initial subject of my internship (by François Fages) here.

Principal modifications of BIOCHAM

The new syntax and how to use it.

The new syntax I introduced allows declaration of molecules' locations and remains compatible with the old one.
First you declare the possible location(s) for the molecules (but all the molecules must not necessary be localised) . Then in the reaction rules, molecules can appear in the localised or in the non-localised form. You can also declare volume ratio between the locations. This is important for numerical simulations when transports occur between two locations.

All the new BIOCHAM commands I have defined is available here. I also add all the changes I have made in my version of the BIOCHAM reference manual. In this version, all my modifications are written in blue.

Molecules Localisation.

In the new syntax, a BIOCHAM object can be a localised molecule:

object = molecule | molecule::name | ...

Now a molecule localised in several compartments can be considered differently according to its location. For instance, Mdm2::cytoplasm is a valid BIOCHAM notation for the Mdm2 protein localised in the cytoplasm.

The set of all possible localisations for a given molecule can be declared by associating to the molecule the set of its possible locations or its single possible localisation.
For example, the location declarations localise p53::nucleus.and localise p53::[nucleus]. denote the same thing: the p53 molecule can be localised in the nucleus (only) . On the other hand localise Mdm2::[cytoplasm,nucleus]. denotes that the Mdm2 molecule can be localised in the cytoplasm or in the nucleus, but in no other specified location.
This molecule could appear in the rules in the localised forms Mdm2::cytoplasm or Mdm2::nucleus or in the non-localised form Mdm2

To use a localised form of a molecule in a rule, this form must have been previously defined with the command localise. Once the locations of a molecule have been specified, you can add new possible locations with the command extend_localise , which is used with the same nomenclature as localise .
All the localised molecules can be listed with the command list_locations and cleared with the command clear_rules.

The location notation is compatible with the old command list_molecules. For instance, if you declare :

localise A::nucleus.
localise B::[nucleus,cytoplasm].
localise C::cytoplasm.
Then you can ask which molecules have been localised in the cytoplasm, or where the B molecule has been declared with the BIOCHAM commands:

Volume ratio.

When several cell compartments are considered, these compartments may have different volumes. In this case, transport rules occur between compartments whose volume differ from each other and volume ratio must be taken into account in all kinetic models.

For instance, the volume ratio between nucleus and cytoplasm (nuclear volume/cytoplasmic volume) is currently above 1/15. It can be defined with the command volume_ratio (15,nucleus),(1,cytoplasm). The volume ratio will automatically be considered during numerical simulations for every rule recognized as transport rules, i.e. every reaction whose at least one reactant and one product are localised in different compartments. If no volume ratio are defined, the default ratio is set to unity (that means all the compartments have the same volume) and a warning is signaled.
All the volume ratio can be listed with the command list_volumes.
If no volume ratio have been declared, the default ratio is 1 for every previously defined location.

These ratio are used for transport rules' numerical simulations.

Transport rules.

In BIOCHAM, rules are written with stoechiometric coefficients, for instance like in the reaction:

A + B => A-B
Here one A molecule and one B molecule give one complex A-B. During numerical simulations, concentrations and molecules numbers can be confused if the volume remains stable during the reaction (Molecules Number = Concentration x Volume).
Of course it isn't a priori the case for transport rules since compartments may be different between right and left member of the rule. During a transport, the compartment volume may change. That's why volume ratio are so important in numerical simulation.
In my BIOCHAM space version, volume ratio are used only for transport rules. Those rules are recognised, following the assumptions:
  • A transport rule is a rule whose at least one reactant of both sides is located,
  • The right and the left side locations are different.

    For instance, those rules are recognised as transport rules:
    A::nucleus => A::cytoplasm
    A~{p}::nucleus => A::cytoplasm (*)
    C::nucleus + B::nucleus => A::cytoplasm
    But these reactions are not:
    A::nucleus => B::nucleus
    C + B::nucleus => A
    C + B::nucleus => B
    (*): Note: A~{p} is the phosphorylated form of the molecule A.

    To avoid confusion every recognised transport is signaled to the user. And when no volume ratio have been defined, the simulation occurs with a default ratio of 1, and a warning is signaled.

    The p53/Mdm2 Network is a good example of the utilisation of volume ratio and of the transport rules recognition.

    Two examples to illustrate the importance of space modeling in BIOCHAM

    Initially, I was supposed to work on the example of the RTK/MAPK signaling cascade. Indeed in the cell the molecules of this pathway are localised: the ligant (EGF for instance) is extracellular, it binds a membrane receptor (RTK), receptors activation activates internal signal: some proteins of the internal membrane side are phosporylated. They activate successively cytoplasmic enzyms. The last enzym of the cascade (ERK) enters the nucleus and regulates gene expression.
    Space modeling becomes useful when the same molecule can be localised in several locations and when this molecule plays a crucial role. In this pathway only the ERK location changes and informations about the ERK transport from the cytoplasm to the nucleus didn't seem to be well documented (Is the transport reversible? What is the kinetic?...). That's why I sought an other example, in witch the transport phenomena would have a central role.

    A Model of the p53/Mdm2 Network

    The p53/Mdm2 Network is a fluctuating system involved in cell apoptosis. In this network, a protein (Mdm2) is actively transported between the cytoplasm and the nucleus.
    P53 is a transcription activator activated under stress conditions in the mammalian cells. In the presence of damaged DNA, the p53 concentration enhances and may stop the cell cycle and lead to apoptosis (programmed death of the cell).
    In normal conditions, p53 is degraded after its ubiquitination by Mdm2 protein, which plays a central role in the p53 regulation and whose transcription is also induced by p53. This negative feedback makes possible to keep low p53 concentration in the normal cells. In cancers, p53 mutation is often observed.
    Experimental studies have shown that in presence of damaged DNA (after irradiation for instance), the p53 and Mdm2 concentrations oscillate a few times and DNA is progressively repair. One mathematical model [2] explains oscillations in the p53/Mdm2 network by the existence of positive and negative feedbacks. It is this model that I have simulated with BIOCHAM. (See original paper for reaction schema)

    Simulation graph
    BIOCHAM simulation of irradiation experiment: p53/Mdm2 oscillations and DNA repair.
    IR: irradiation; Mdm-Mdm::n : nuclear Mdm2 protein;
    Mdm-Mdm::c : cytoplasmic Mdm2 protein; ;
    Mdm-Mdm~{p}::c : phosphorylated cytoplasmic Mdm2 protein
    p53, p53-u, p53-u-u: non-, mono- and bi-ubiquitinated p53 protein.

    When DNA is damaged (blue line), two oscillations of p53 and Mdm2 concentrations occur.
    The damaged DNA activates Mdm2 degradation, what raises the inhibition on p53 ubiquitination: p53 is less degraded. In one hand its presence in the cell leads to the repair of DNA, in the other hand p53 catalyses the Mdm2 production in the cytoplasm, and Mdm2 is transported to the nucleus, and activates p53 ubiquitination and degradation.

    Mdm2 transport: importance of the volume ratio

    The Mdm2 transport through the nuclear envelope plays the key role of the networks regulation. Once phosphorylated the cytoplasmic form of Mdm2 can be actively transported to the nucleus (and dephosphorylated). In BIOCHAM localised syntax the reaction rule (without kinetics coefficients) is:

    localise Mdm-Mdm::[cytoplasm,nucleus].
    localise Mdm-Mdm~{p}::cytoplasm.

    add_rule(Mdm-Mdm::nucleus <=> Mdm-Mdm~{p}::cytoplasm).

    In the localised version, this reaction rule is recognised as a transport rule. If no volume ratio is declared (a warning is signaled), the numerical simulation of this reaction is wrong, because even if one nuclear Mdm2 molecule gives one cytoplasmic (and reciprocally), the ratio don't stay true for the concentrations because the cytoplasmic volume is about 15 times bigger then the nuclear volume.
    In order to be taken into account the volume ratio must be previously defined with the command:
    volume_ratio (15,nucleus),(1,cytoplasm).

    The Delta-Notch signaling: a funny example about a possible use of BIOCHAM in the future

    Delta and Notch are two proteins involved in the cellular differentiation. Their synthesis depends on their presence in the neighbour cells. The hybrid model proposed by [3] explains how, with lateral inhibitions, you can obtain the salt-and-paper pattern, observed in biological experiments:

  • At the steady state, all cells have either the Delta phenotype, or the Notch.
  • No Delta cell has a Delta neighbour, and a Notch cell can't only have Notch neighbours.

    In the model, all initial concentrations are taken randomly to a normal distribution. Then the Delta concentration of each cells varies according to the Notch concentration of the cell, and the Notch concentration of each cells varies according to the Delta concentrations of its neighbours. What leads to a steady state.
    The simulation has been tested for an initial state close to a normal distribution and with several configurations (50 cells loop and 2 neighbours per cell, 6 x 6 cells pattern and 4 neighbours per cell).

    Here is a diagram of the Delta/Notch differentiation for a 6 x 6 cells repetitive pattern, with a model with 4 neighbours per cell (the steady state is reached after ten time units).

    Delta/Notch Differentiation in a repetitive 6 x 6 cells array.
    Cell phenotypes: Delta: blue, Notch : yellow.

    A possible development for the future

    To use location relations to simplify rule writing

    It would be interesting to declare location relations like neighbourhood or non-accessibility in order to simplify rule writing. When locations are declared, it would be easy to declared them with neighbourhood classes (set of 'neighbour' locations). For instance, the cytoplasm and the nucleus may be considered as 'neighbours locations' in the Mdm2 Network example, it's of course also the case for neighbours cells in the Delta/Notch signaling.

    With the Delta/Notch example, it would allow to generalize rule writing since the Notch synthesis is conditioned only by the neighbour cells.

    To guess the locations and their relations

    Explicit localisation

    In the spacial version I have implemented, the molecule localisation must be explicit.
    Initially, I considered the possibility to extend the locations for the rules. I. e. that if the molecules are declared in one location:

    localise A::L1.
    localise B::L2.
    localise C::L3.
    and if you consider the rule where A + B => C then BIOCHAM would automaticaly localise the rule in A::L1 + B::L2 => C::L3 . But, how localised the rule if one or several molecules are localised in more than one location? That's why I chose to impose to the user to write explicitly the localised form in the rules.
    Even if a rule occurs in one location, you have to write all molecules in their localised form.

    Generalized localisation

    If a reaction occurs in several locations, you also have to write the rules for each location. It would be useful to develop a version where the locations may be implicit in such a case. For instance if the rule A + B => C occurs in location L1 and in location L2, in the current version, you have to write:

    localise A::[L1,L2].
    localise B::[L1,L2].
    localise C::[L1,L2].

    add_rule(A::L1 + B::L1 => C::L1).
    add_rule(A::L2 + B::L2 => C::L2).
    It would be a possibility to write these two commands in one, like:
    add_rule(A + B => C,L1,L2).

    In the Delta/Notch signaling for instance, the rules of Delta and Notch synthesis are the same for each cell. With the current BIOCHAM space syntax, you have to write explicitly the rules for all the cells you model. With a lighter syntax, it would be more simple for users, especially if neighbourhood classes are declared.

    An other possible extention, would be to guess locations and their relations from the rules, according the assumption that if several locations appear in one rule then these locations are 'neighbours'.

    To model intracellular reactions, the location relation may seem superfluous if only a few locations are considered (extracellular domain, plasmic membrane, cytoplasm, nuclear envelop, nucleus). But to model intercellular reactions, a lot of locations (cells) may be considered and such an extention would be useful.


    I would like to thank François Fage, for giving me the oportunity to da this internship and for welcoming me in his team.

    Thanks to all the Contraintes team for its cordial welcome, especially Nathalie and Laurence for supporting me in their office and for their precious help.

    I would also like to thank my friends for helping me and contributing to make this internship pleasant, and in general everyone I had the occasion to meet, at the coffee machine for instance.


    [1] The Biochemical Abstract Machine BIOCHAM by Nathalie Chabrier, François Fages and Sylvain Soliman.
    Proceedings of the Second International Workshop on Computational Methods in Systems Biology CMSB'04, Springer-Verlag, Paris, France, May 2004.
    Available as ps.gz.

    [2] Steady States and Oscillations in the p53/Mdm2 Network by Andrea Ciliberto, Béla Novak, John J. Tyson.
    Cell Cycle 4:3, 488-493; March 2005.

    [3] Lateral inhibition through Delta-Notch signaling: A piecewise affine hybrid model by R. Ghosh and C. Tomlin.
    Hybrid Systems: Computation ann Control, HSCC 2001.Available as pdf.

    [4] Using hybrid concurrent constraint programming to model dynamic biological systems by A. Bockmayer and A. Courtois. Techreport, Université Henry Poincaré, LORIA, Vandoeuvre-les-Nancy, France, 2002. Available as ps.

    Sakina Ayata
    Last modified: Fri Jun 3 16:51:56 CEST 2005