Luca Cardelli Microsoft Research Cambridge UK Cambridge 2005-07-06 Anatomy Department www.luca.demon.co.uk/BioComputing.htm

1
Luca CardelliMicrosoft ResearchCambridge
UKCambridge 2005-07-06Anatomy
Departmentwww.luca.demon.co.uk/BioComputing.htm
Luca CardelliMicrosoft ResearchCambridge
UKSerious Talk, Cambridge 2005-03-23www.luca.d
emon.co.uk
Biological Systems as Reactive Systems
2
50 Years of Molecular Cell Biology

• Genes are made of DNA
• Store digital information as sequences of 4
  different nucleotides
• Direct protein assembly through RNA and the
  Genetic Code
• Proteins (gt10000) are made of amino acids
• Process signals
• Activate genes
• Move materials
• Catalyze reactions to produce substances
• Control energy production and consumption
• Bootstrapping still a mystery
• DNA, RNA, proteins, membranes are today
  interdependent. Not clear who came first
• Separation of tasks happened a long time ago
• Not understood, not essential

3
Towards Systems Biology

• Biologists now understand many of the cellular
  components
• A whole team of biologists will typically study a
  single protein for years
• Reductionism understand the components in order
  to understand the system
• But this has not led to understand how the
  system works
• Behavior comes from complex patterns of
  interactions between components
• Predictive biology and pharmacology still rare
• Synthetic biology still unreliable
• New approach try to understand the system
• Experimentally massive data gathering and data
  mining (e.g. Genome projects)
• Conceptually modeling and analyzing networks
  (i.e. interactions) of components
• What kind of a system?
• Just beyond the basic chemistry of energy and
  materials processing
• Built right out of digital information (DNA)
• Based on information processing for both survival
  and evolution
• Highly concurrent

4
Storing Processes

• Today we represent, store, search, and analyze
• Gene sequence data
• Protein structure data
• Metabolic network data
• Signaling pathway data
• How can we represent, store, and analyze
  biological processes?
• Scalable, precise, dynamic, highly structured,
  maintainable representations for systems biology.
• Not just huge lists of chemical reactions or
  differential equations.
• In computing
• There are well-established scalable
  representations of dynamic reactive processes.
• They look more or less like little,
  mathematically based, programming languages.

Cellular Abstractions Cells as
Computation RegevShapiro NATURE vol 419,
2002-09-26, 343
5
Structural Architecture
Nuclear membrane
EukaryoticCell (10100 trillion in human body)
Mitochondria
Membranes everywhere
Golgi
Vesicles
E.R.
Plasma membrane (lt10 of all membranes)
H.Lodish et al. Molecular Cell Biology fourth
edition p.1
6
Abstract Machines of Systems Biology
Regulation
The hardware (biochemistry) is fairly well
understood.But what is the software that runs
on these machines?
GeneMachine
Gene Regulatory Networks
strings
Nucleotides
Functional ArchitectureDiverse - chemical
toolkits- instruction sets- programming
models- notations
Makes proteins,where/when/howmuch
Holds genome(s),confines regulators
Directs membrane construction and protein
embedding
Signals conditions and events
Model Integration Different time and space scales
Holds receptors, actuators hosts reactions
Transport Networks
ProteinMachine
Machine
Biochemical Networks
Membrane
Implements fusion, fission
Aminoacids
records
Phospholipids
hierarchical multisets
Metabolism, PropulsionSignal ProcessingMolecular
Transport
ConfinementStorageBulk Transport
trees
7
Reactive Systems
Probably complexity is in not in any fundamental
sense rarer in continuous systems than in
discrete ones. But the point is that discrete
systems can typically be investigated in a much
more direct way than continous ones. Stephen
Wolfram, A New Kind of Science, p.167

• Modeling biological systems
• Not as continuous systems (often highly
  nonlinear)
• But as discrete reactive systems abstract
  machines where
• States represent situations
• Event-driven transitions between states represent
  dynamics
• The adequacy of describing (discrete) complex
  systems as reactive systems has been argued
  convincingly David Harel
• Many biological systems exhibit features of
  reactive systems
• Discrete transitions between states
• Deep layering of abstractions (steps at
  multiple levels)
• Complexity from combinatorial interaction of
  simple components
• High degree of concurrency and nondeterminism
• Emergent behavior not obvious from part list
• Still, needs quantitative semantics
• Stochastic, hybrid, etc. to talk about rates (and
  geometry).

8
Methods

• Model Construction (writing things down
  precisely)
• Formalizing the notations used in systems
  biology.
• Formulating description languages.
• Studying their kinetics (semantics).
• Model Validation (using models for postdiction
  and prediction)
• Simulation from compositional descriptions
• Stochastic quantitative concurrent semantics.
• Hybrid discrete transitions between continuously
  evolving states.
• Program Analysis
• Control flow analysis
• Causality analysis
• Modelchecking
• Standard, Quantitative, Probabilistic

9
Chemistry vs. p-calculus
A process calculus (chemistry)
A different process calculus (p)
r A B ?k1 C Ds C D ?k2 A B
A
B
Does A become C or D?
A
B
!rk1
?rk1
?sk2
!sk2
Reactionoriented
1 line per reaction
rk1
C
D
Interactionoriented
Interactionoriented
1 line per component
C
D
A !rk1 CC ?sk2 AB ?rk1 DD
!sk2 B
A becomes C not D!
sk2
A Petri-Net-like representation. Precise and
dynamic, but not modular, scalable, or
maintainable.
A compositional graphical representation
(precise, dynamic and modular) and the
corresponding calculus.
10
Biochemistry vs. p-calculus
A, B, C
ABC 69 ApBC ABC 69 ABpC ABC 69 ABCp ApBC 69
ApBpC ApBC 69 ApBCp ABpC 69 ApBpC ABpC 69
ABpCp ABCp 69 ApBCp ABCp 69 ABpCp ApBpC 69
ApBpCp ApBCp 69 ApBpCp ABpCp 69 ApBpCp
A 69 Ap B 69 Bp C 69 Cp
domains
A 69 Ap B 69 Bp C 69 Cp
ABC ApBC ABpC ABCp ApBpC ApBCp ABpCp ApBpCp
2n
domainreactions
domainreactions
2n(2n-1)
ABC
A ?knAp Ap ?phA B ?knBp Bp ?phB C
?knCp Cp ?phC A B C
reactions(twice number of edges in n-dim
hypercube)
species
complex
processes
Stoichiometric Matrix (species x reaction)System
of Differential Equations
The matrix is very sparse, so the corresponding
ODE system is not dense. But it still has 2n
equations, one per species, plus conservation
equations (ABCApBCconstant, etc.).
System description is linear in the number of
basic components.
System description is exponential in the number
of basic components.
2n x 2n(2n-1)
(Its run-time behavior or analysis potentially
blows-up just as in the previous case, but its
description does not.)
(twice number of edges in n-dim hypercube, which
is given by Sloanes Sequence A001787
0,1,4,12,32,80,192,448,1024,2304,5120,11264,24576,
53248, 114688,245760,524288,1114112,2359296,498073
6,10485760, 22020096,46137344,96468992,201326592,4
19430400,872415232, 1811939328,3758096384,77846282
24 )
When the local domain reactions are not
independent, we can use lateral communication so
that each component is aware of the relevant
others.
11
Why p-calculus, in particular

• Well studied, compact, precise, and general
• A programming language first, a mathematical
  model second
• Syntax (configurations) P 0 PP PP
  ?n(n).P !n(n).P (nn)P P
• Semantics (reactions) (P !n(m).P) (Q
  ?n(m).Q) P Qmm
• Binary interactions
• I.e., collisions
• Reactive and compositional
• Each subsystem is a separate (composition of)
  automata interacting with the environment (more
  automata)
• Dynamic network evolution and species evolution
• Each subsystem can create fresh connections or
  spawn new subsystems
• Compact description of combinatorics (like any
  programming language)
• (Bit1 Bit2 Bitn) size-n description,
  where Bit is a 2-state subsystem
• E.g. a protein with with 2n phosphorilation
  configurations (i.e. different chemical
  species)
• Complexation/polymerization
• The most characteristic feature of p-calculus
  (fresh names) models sticking

12
Basic Modeling Guidelines

• Regev-Shapiro Molecules as Processes
• They chose p-calculus and adapted it with
  stochastic features
• To match the stochastic aspects of (bio)chemistry
• Many probabilistic process calculi predate them,
  but only Hillston (CSP) and Priami (p) had
  already studied stochastic calculi.

Molecule Process
Interaction capability Channel
Interaction Communication
Modification (of chemical components) State change(state-transition systems)

• Cellular Abstractions Cells as Computation
• RegevShapiro NATURE vol 419, 2002-09-26, 343

13
p-calculus Executive Summary

• Its for
• The modular description of concurrent,
  nondeterministic systems
• Study of such systems based on their descriptions
• Its got
• Processes
• Channels
• A minimalistic syntax (its a language and also a
  model)
• You can
• Fork new processes
• Create new channels
• Do I/O over channels (synchronous and
  asynchronous)including passing channels over
  channels
• Make nondeterministic choices
• Define processes recursively
• Thats it.
• Except for extensive model theory and metatheory.
• Cannot pass processes over channels (simulated
  by passing channels to them)

14
p-calculus (a Process Algebra)

• Processes P,Q, - components of a system
• Channels a,b, - interactions between
  components
• 0 the process that does nothing
• !a(b) P the process that outputs b on channel a
  (and then does P)
• ?a(x) P the process that inputs b on channel a
  (and then does Px)
• P Q the process made of subprocesses P and Q
  running concurrently
• P Q the process that behaves like either P or Q
  nondeterministically
• P the process that behaves like unboundedly
  many copies of P
• gt recursive processes
• gt unbounded number and species of processes
• new x P the process that creates a new channel x
  (and then does Px)
• gt private interactions
• gt unbounded number and species of interactions

15
p-calculus (a Process Algebra)

• Dynamics
• (!a(b) P) P (?a(x) Qx) Q
  ? P Qb
• Compositional descriptions
• Describe how the individual components behave
• i.e. how they interact with any environment they
  may be placed in
• Build systems by combining components
• each components is part of the environment for
  the other components
• Behavior (and its analysis) arises from the
  combinatorics of interactions
• state space can be arbitrarily larger than its
  compositional description
• For concurrent, nondeterministic, unbounded-state
  systems
• Dynamic creation of new channels (e.g. binding
  sites)
• Dynamic creation of new processes (e.g. proteins)

16
p-calculus
Syntax
Chemical Mixing
Reactions
17
Stochastic p-calculus

• A quantitative variant of p-calculus
• Channels have stochastic firing rates with
  exponential distribution.
• Nondeterministic choice becomes stochastic race.
• Cuts down to CTMCs (Continuous Time Markov
  Chains) in the finite case (not always). Then,
  standard analytical tools are applicable.
• Can be given friendly automata-like scalable
  graphical syntax (Andrew Phillips et al.).
• Is directly executable (e.g. via the Gillespie
  algorithm from physical chemistry).
• Is analyzable (large body of literature, at least
  in the non-stochastic case).

A.Phillips, L.Cardelli. BioConcur04.
18
Importance of Stochastic Effects

• A deterministic system
• May get stuck in a fixpoint.
• And hence never oscillate.
• A similar stochastic system
• May be thrown off the fixpoint by stochastic
  noise, entering a long orbit that will later
  bring it back to the fixpoint.
• And hence oscillate.

Mechanisms of noise-resistance in genetic
oscillators Jose M. G. Vilar, Hao Yuan Kueh,
Naama Barkai, Stanislas Leibler PNAS April 30,
2002 vol. 99 no. 9 p.5991
19
Gene Networks
20
The Gene Machine
The Central Dogma of Molecular Biology
regulation
transcription
translation
interaction
folding
21
The Gene Machine Instruction Set
cf. Hybrid Petri Nets Matsuno, Doi, Nagasaki,
Miyano
Positive Regulation
Transcription
Negative Regulation
Input
Output
Coding region
Gene(Stretch of DNA)
External Choice The phage lambda switch
Regulatory region
Regulation of a gene (positive and negative)
influences transcription. The regulatory region
has precise DNA sequences, but not meant for
coding proteins meant for binding
regulators. Transcription produces molecules (RNA
or, through RNA, proteins) that bind to
regulatory region of other genes (or that are
end-products).
Human (and mammalian) Genome Size3Gbp (Giga base
pairs) 750MB _at_ 4bp/Byte (CD) Non-repetitive
1Gbp 250MB In genes 320Mbp 80MB Coding
160Mbp 40MB Protein-coding genes
30,000-40,000 M.Genitalium (smallest true
organism) 580,073bp 145KB (eBook)E.Coli
(bacteria) 4Mbp 1MB (floppy)Yeast (eukarya)
12Mbp 3MB (MP3 song)Wheat 17Gbp 4.25GB (DVD)
22
Gene Composition
Is a shorthand for
a
b
Under the assumptions Kim Tidor1) The
solution is well-stirred (no spatial
dependence on concentrations or rates).2) There
is no regulation cross-talk.3) Control of
expression is at transcription level only
(no RNA-RNA or RNA-protein effects)4)
Transcriptions and translation rates
monotonically affect mRNA and protein
concentrations resp.
Ex Bistable Switch
a
b
a
b
Ex Oscillator
Expressed
c
c
c
Repressed
Expressing
a
b
a
b
a
b
23
Gene Regulatory Networks
http//strc.herts.ac.uk/bio/maria/NetBuilder/
NetBuilder
24
(The Classical ODE Approach)
Chen, He, Church
I.e. to model an operating system, write a set
of differential equations relating the
concentrations in memory of data structures and
stack frames over time. (Duh!)
n number of genesr mRNA concentrations (n-dim
vector)p protein concentrations (n-dim
vector)f (p) transcription functions (n-dim
vector polynomials on p)
L r - U r
25
Nullary Gate
spontaneous (constitutive) output
b
no input
null
interaction site of output protein
null(b) _at_ te (tr(b) null(b))
(recursive, parametric) process definition
and repeat
output protein (transcripion factor), spawn out
stochastic delay (t) with rate e of constitutive
transcription
A stochastic rate r is always associated with
each channel ar (at channel creation time) and
delay tr, but is often omitted when unambiguous.
26
Production and Degradation
Degradation is extremely important and often
deliberate it changes unbounded growth into
(roughly) stable signals.
and repeat
transcripton factor
degradation
tr(p) _at_ (!pr tr(p)) td
degradation rate d
(output, !) interaction with rate r (input, ?, is
on the target gene)
interaction site of transcription factor
stochastic choice (race between r and d)
A transcription factor is a process (not a
message or a channel) it has behavior such as
interaction on p and degradation.
combined effect of production and degradation
(without any interaction on b)
null(b)
e0.1, d0.001
b
product
interaction offers on b ( number of tr processes)
b
null(b) _at_ te (tr(b) null(b))
null
time
27
Unary Pos Gate
output (stimulated or constitutive)
input (excitatory)
transcripton delay with rate h
pos(a,b) _at_ ?ar th (tr(b) pos(a,b))
te (tr(b) pos(a,b))
(input, ?) interaction with rate r
race between r and e
or constitutive transcription to always get
things started
output protein
parallel, not sequence, to handle self-loops
without deadlock
unlimited amount of
r1.0, e0.01, h0.1, d0.001
b
Stimulated
tr(ar) pos(ar,b)
pos(a,b)
Constitutive
28
Unary Neg Gate
output (constitutive when not inhibited)
input (inhibitory)
inhibition delay with rate h
neg(a,b) _at_ ?ar th neg(a,b) te (tr(b)
neg(a,b))
(input, ?) interaction with rate r
or constitutive transcription to always get
things started
race between r and e
r1.0, e0.1, h0.01, d0.001
b
Constitutive
neg(ar,b)
tr(ar) neg(ar,b)
Inhibited
29
Signal Amplification
pos(a,b) _at_ ?ar th (tr(b) pos(a,b))
te (tr(b) pos(a,b))
E.g. 1 a that interacts twice before decay can
produces 2 b that each interact twice before
decay, which produce 4 c
pos(a,b) pos(b,c)
a
c
b
pos
pos
tr(p) _at_ (!pr tr(p)) td
30
Signal Normalization
neg(a,b) _at_ ?ar th neg(a,b) te (tr(b)
neg(a,b))
neg(a,b) neg(b,c)
a
c
b
neg
neg
tr(p) _at_ (!pr tr(p)) td
r1.0, e0.1, h0.01, d0.001
a non-zero input level, a, whether weak or
strong, is renormalized to a standard level, c.
b
c
a
30tr(a) neg(a,b) neg(b,c)
31
Self Feedback Circuits
pos(a,a)
neg(a,a)
a
a
neg
pos
neg(a,b) _at_ ?ar th neg(a,b) te (tr(b)
neg(a,b))
pos(a,b) _at_ ?ar (tr(b) pos(a,b)) te
(tr(b) pos(a,b))
tr(p) _at_ (!pr tr(p)) td
tr(p) _at_ (!pr tr(p)) td
(Can overwhelm degradation, depending on
parameters)
high, to raise the signal
r1.0, e10.0, h1.0, d0.005
a
neg(a,a)
32
Two-gate Feedback Circuits
pos(b,a) neg(a,b)
neg(b,a) neg(a,b)
Bistable
Monostable
For some degradation rates is quite stable
r1.0, e0.1, h0.01, d0.001
a
b
a
b
neg(b,a) neg(a,b)
But with a small change in degradation, it goes
wild
e0.1, h0.01, d0.001
r1.0, e0.1, h0.01, d0.0001
a
5 runs with r(a)0.1, r(b)1.0 shows that circuit
is now biased towards expressing b
b
b
pos(b,a) neg(a,b)
33
Repressilator
neg(a,b) _at_ ?ar th neg(a,b) te (tr(b)
neg(a,b))
neg(a,b) neg(b,c) neg(c,a)
Same circuit, three different degradation models
by chaning the tr component
interact once and die otherwise stick around
interact once and die otherwise decay
tr(p) _at_ !pr
tr(p) _at_ !pr td
r1.0, e0.1, h0.04
r1.0, e0.1, h0.04, d0.0001
a b c
a b c
interact many times and decay
tr(p) _at_ (!pr tr(p)) td
r1.0, e0.1, h0.001, d0.001
a b c
Subtle at any point one gate is inhibited and
the other two can fire constitutively. If one of
them fires first, nothing really changes, but if
the other one fires first, then the cycle
progresses.
34
System Properties Oscillation Parameters
The constitutive rate e (together with the
degradation rate) determines oscillation
amplitude, while the inhibition rate h determines
oscillation frequency.
We can view the interaction rate r as a measure
of the volume (or temperature) of the solution
that is, of how often transcription factors bump
into gates. Oscillation frequency and amplitude
remain unaffected in a large range of variation
of r.
35
Repressilator in SPiM
val dk 0.001 ( Decay rate ) val eta
0.001 ( Inhibition rate ) val cst 0.1 (
Constitutive rate ) let tr(pchan()) do !p
tr(p) or delay_at_dk let neg(achan(), bchan())
do ?a delay_at_eta neg(a,b) or delay_at_cst
(tr(b) neg(a,b)) ( The circuit ) val bnd
1.0 ( Protein binding rate ) new a_at_bnd
chan() new b_at_bnd chan() new c_at_bnd
chan() run (neg(c,a) neg(a,b) neg(b,c))
36
System Properties Fixpoints
A sequence of neg gates behaves as expected, with
alternating signals, (less Booleanly depending
on attenuation).
Now add a self-loop at the head. Not a Boolean
circuit!No more alternations, because each
gate is at its fixpoint.
unstable
all low!
37
Repressilator ODE Model and Simulation
Bruce E Shapiro Cellerator
38
Guet et al. D038/lac-
Combinatorial Synthesis of Genetic Networks,
Guet, Elowitz, Hsing, Leibler, 1996, Science, May
2002, 1466-1470.
IPTG
aTc
The output of some circuits did not seem to make
any sense
D038/lac-
lcI
GFP
TetR
LacI
PT
PL2
PT
Pl-
tet
lac
cI
gfp
neg(TetR,TetR) neg(TetR,LacI) neg(LacI,lcI)
neg(lcI,GFP)
r1.0, e0.1, h1.0, d0.001
aTc -IPTG -GFP -
aTc -IPTG GFP -
aTc IPTG GFP -
aTc IPTG -GFP
GFP!
A Compositional Approach to the Stochastic
Dynamics of Gene Networks, Ralf Blossey, Luca
Cardelli, Andrew Phillips, TCSB, Springer, to
appear.
39
Guet et al. D038/lac-
Combinatorial Synthesis of Genetic Networks,
Guet, Elowitz, Hsing, Leibler, 1996, Science, May
2002, 1466-1470.
r1.0, e0.1, h1.0, d0.001
aTc -IPTG -GFP -
aTc -IPTG GFP -
We can model an inducer like aTc as something
that competes for the transcription factor.
IPTG de-represses the lac operon, by binding to
the lac repressor (the lac I gene product),
preventing it from binding to the operator.
aTc IPTG -GFP
GFP!
aTc IPTG GFP -
40
Guet et al.
Combinatorial Synthesis of Genetic Networks,
Guet, Elowitz, Hsing, Leibler, 1996, Science, May
2002, 1466-1470.
They engineered in E.Coli all genetic circuits
with four single-input gates such as this one
We can model an inducer like aTc as something
that competes for the transcription factor.
Then they measured the GFP output (a fluorescent
protein) in presence or absence of each of two
inhibitors (aTc and IPTG).
The output of some circuits did not seem to make
any sense
IPTG de-represses the lac operon, by binding to
the lac repressor (the lac I gene product),
preventing it from binding to the operator.
Here 1 means high brightness and 0 means
low brightness on a population of bacteria
after some time. (I.e. integrated in space and
time.)
41
Further Building Blocks
42
D038/lac-
Naïve Boolean analysis would suggest GFP0.5
(oscillation) because of self-loop.
GFP0 there is consistent only with (somehow) the
head loop setting TetRLacI0. But in that case,
aTc should have no effect (it can only subtract
from those signals) but instead it sets GFP1.
Hence we need to understand better the dynamics
of this network.
43
Simulation results for D038/lac-
We can model an inducer like aTc as something
that competes for the transcription factor.
IPTG de-represses the lac operon, by binding to
the lac repressor (the lac I gene product),
preventing it from binding to the operator.
44
D016/lac-
How can aTc affect the result??
One theory aTc prevents the self-inhibition of
tet, so that a very large quantity of TetR is
produced. That then overloads the overall
degradation machinery of the cell, affecting the
rest of the circuit.
Even so, how can GFP be high here?
Even the fixpoint explanation fails here, unless
we assume that the lac gate is operating in its
instability region.
45
Simulation results for D016/lac-
A
B
aTc 1 (d 0.00001), IPTG 0
GFP
The fixpoint effect, in instability region,
explains this GFP high because wildly
oscillating.
The fixpoint effect, in instability region,
explains this GFP high because wildly
oscillating.
C
D
aTc 0 (d 0.001), IPTG 1
aTc 1 (d 0.00001), IPTG 1
Overloading of degradation machinery, induced by
aTc, can reinstate the fixpoint regime.
Overloading of degradation machinery, induced by
aTc, can reinstate the fixpoint regime.
E
r 1.0e 0.1h 0.01
d 0.005 aTc 0, IPTG 0
46
What was the point?

• Deliberately pick a controversial/unsettled
  example to test the methodology.
• Show that we can easily play with the model and
  run simulations.
• Get a feeling for the kind of subtle effects that
  may play a role.
• Get a feeling for kind of analysis that is
  required to understand the behavior of these
  systems.
• In the end, we are never understanding
  anything we are just building theories/models
  that support of contradict experiments (and that
  suggest further experiments).

47
Protein Networks
Luca CardelliMicrosoft ResearchCambridge
UK2005-07www.luca.demon.co.uk
48
MAPK Cascade - HuangFerrell
Ultrasensitivity in the mitogen-activated protein
cascade, Chi-Ying F. Huang and James E. Ferrell,
Jr., 1996, Proc. Natl. Acad. Sci. USA, 93,
10078-10083.
10 chemical reactions
Reservoirs
Reservoirs
Reservoirs
Back Enzymes
Back Enzymes
49
As 18 Ordinary Differential EquationsPlus 7
conservation equations
One equation for each species (8) and complex
(10), but not for constant concentration enzymes
(4)
in exactly one state
Each molecule
50
The Circuit
51
Enzymatic Reactions
Reaction View
intermediatecomplex
º
E
(c,d,e)
S
P
Interaction View
private bindings between one S and one E
molecule
bind
S() _at_ new u_at_d new k_at_e !ac(u,k) (!ud S()
!ke P())
unbind
react
bind
unbind
react
ac
ud
ke
E() _at_ ?ac(u,k) (?ud E() ?ke E())
P() _at_
52
MAPK Cascade in SPiM
let KKK() (new u1_at_d1Release new
k1_at_r1React !a1(u1,k1) (do !u1KKK() or
!k1KKKst())) and KKKst() (new
u2_at_d2Release new k2_at_r2React do !a2(u2,k2)
(do !u2KKKst() or !k2KKK()) or
?a3(u3,k3) (do ?u3KKKst() or ?k3KKKst())
or ?a5(u5,k5) (do ?u5KKKst() or
?k5KKKst())) let E1() ?a1(u1,k1) (do
?u1E1() or ?k1E1()) let E2() ?a2(u2,k2)
(do ?u2E2() or ?k2E2()) let KK() (new
u3_at_d3Release new k3_at_r3React !a3(u3,k3)
(do !u3KK() or !k3KK_P())) and KK_P()
(new u4_at_d4Release new k4_at_r4React new
u5_at_d5Release new k5_at_r5React do !a4(u4,k4)
(do !u4KK_P() or !k4KK()) or !a5(u5,k5)
(do !u5KK_P() or !k5KK_PP()))
and KK_PP() (new u6_at_d6Release new
k6_at_r6React do !a6(u6,k6) (do !u6KK_PP()
or !k6KK_P()) or ?a7(u7,k7) (do
?u7KK_PP() or ?k7KK_PP()) or ?a9(u9,k9)
(do ?u9KK_PP() or ?k9KK_PP())) and KKPse()
do ?a4(u4,k4) (do ?u4KKPse() or
?k4KKPse()) or ?a6(u6,k6) (do ?u6KKPse()
or ?k6KKPse()) let K() (new
u7_at_d7Release new k7_at_r7React !a7(u7,k7)
(do !u7K() or !k7K_P())) and K_P()
(new u8_at_d8Release new k8_at_r8React new
u9_at_d9Release new k9_at_r9React do
!a8(u8,k8) (do !u8K_P() or !k8K()) or
!a9(u9,k9) (do !u9K_P() or !k9K_PP())) and
K_PP() (new u10_at_d10Release new
k10_at_r10React !a10(u10,k10) (do
!u10K_PP() or !k10K_P())) and KPse() do
?a8(u8,k8) (do ?u8KPse() or ?k8KPse()) or
?a10(u10,k10) (do ?u10KPse() or ?k10KPse())
1substrate 2substrate 3kinase 5kinase
1enzyme 2enzyme 3substrate 4su
bstrate 5substrate
6substrate 7kinase 9kinase 4phtase 6
phtase 7substrate 8substrate 9substra
te 10substrate 8phtase 10phtase
One process for each component (12) including
enzymes, but not for complexes.
KKKE1 complex
No need for conservation equations implicit in
choice operator in the calculus.
E1KKK complex
53
globals
type Release chan() type React chan() type
Bond chan(Release,React) new a1_at_1.0Bond val
d11.0 val r11.0 new a2_at_1.0Bond val d21.0 val
r21.0 new a3_at_1.0Bond val d31.0 val r31.0 new
a4_at_1.0Bond val d41.0 val r41.0 new a5_at_1.0Bond
val d51.0 val r51.0 new a6_at_1.0Bond val d61.0
val r61.0 new a7_at_1.0Bond val d71.0 val
r71.0 new a8_at_1.0Bond val d81.0 val r81.0 new
a9_at_1.0Bond val d91.0 val r91.0 new
a10_at_1.0Bond val d101.0 val r101.0 run 100
of KKK() run 100 of KK() run 100 of K() run 1
of E2() run 1 of KKPse() run 1 of KPse() run 1
of E1()
ai(ui,ki) release (ui_at_di) and react (ki_at_ri)
channels passed over bond (ai) channel. (No
behavior attached to channels except interaction
rate.)
54
MAPK Cascade Simulation in SPiM
1st stage KKK barely rises2nd stage KK-PP
rises, but is not stable 3rd stage K-PP flips
up to max even anticipating 2nd stage
K
KK
K-PP
KKK
KK-PP
KKK
Rates and concentrations from paper 1xE2 (0.3
nM) 1xKKPase (0.3 nM) 120xKPase (120 nM) 3xKKK (3
nM) 1200xKK (1.2 uM) 1200xK (1.2 uM) dx rx
150, ax 1 (Kmx (dx rx) / ax, Km 300
nM) 1xE1
KK-P
K-P
55
MAPK Cascade Simulation in SPiM
All coefficients 1.0 !!! 100xKKK, 100xKK, 100xK,
13xE2, 13xKKPse, 13xKPse. nxE1 as
indicated (1xE1 is not sufficient to produce an
output)
56
MAPK Cascade Simulation in SPiM
1st stage KKK barely rises2nd stage KK-PP
rises, but is not stable 3rd stage K-PP flips
up to max even anticipating 2nd stage
KKK
KKK
KK
KK-PP
K
K-PP
All coefficients 1.0 !!! 100xKKK, 100xKK, 100xK,
5xE2, 5xKKPse, 5xKPse. Input is 1xE1.
Output is 90xK-PP (ultrasensitivity).
KK-P
K-P
1xE1 injected
57
Gene-Protein Networks
58
Indirect Gene Effects
No combination of standard high-throughput
experiments can reconstruct an a-priori known
gene/protein network Wagner.
One of many bistable switches that cannot be
described by pure gene regulatory networks
Francois Hakim.
59
François Hakim Fig3A
PNAS (101)2, 580-585, 2004
Design of genetic networks with specified
functions by evolution in silico
Fig 3A
Reactionoriented
Free evolution
Fig 14A
60
François Hakim Fig3A, SPiM simulation
Parameters as in paper
3 copies of each gene.
SPiM simulation
SPiM simulation
SPiM simulation
Spontaneous switch at 500(as discussed in
Supporting Text)30xB injected at 300030xA
injected at 4000
Free evolution
Spontaneous switch at 1100100xB injected at
300030xA injected at 4000
Modified for stability dkA 0.02, dkB 0.02
SPiM simulation
SPiM simulation
SPiM simulation
61
François Hakim Fig3Ast8
Circuit of Fig 3A with parameters from
SupportingText Fig 8, plotted in Fig 13A
SPiM simulation
SPiM simulation
SPiM simulation
SPiM simulation
200xA injected at 2500500xB injected at
5000200xA injected at 7500
200xB injected at 0600xA injected at 2500600xB
injected at 7500
62
François Hakim 3A in SPiM
( Francois and Hakim circuit 3A ) val pntAunb
0.42 val geneACst 0.20 val geneBCst
0.37 val geneBInh 0.027 val bA 0.19 val AB
0.72 val dkA 0.0085 val dkB 0.034 val dkAB
0.53
let ptnA() (new unb_at_pntAunb do
delay_at_dkA or !AB or !bA(unb)(?unb ptnA())) let
ptnB() do delay_at_dkB or ?ABcpxAB() let
cpxAB() delay_at_dkAB let geneA()
delay_at_geneACst (ptnA() geneA()) let
geneBfree() do delay_at_geneBCst (ptnB()
geneBfree()) or ?bA(unb) geneBbound(unb) and
geneBbound(unbch()) do delay_at_geneBInh
(ptnB() geneBbound(unb)) or !unb
geneBfree() run (geneA() geneBfree())
Interactionoriented
63
Towards the Million-Line Model
64
From Chemical Reactions to ODEs
r1 AB ?k1 CC r2 AC ?k2 D r3 C ?k3 EF r4 F
?k4 B
Write the coefficients by columns
Stoichiometric Matrix
N r1 r2 r3 r4
A -1 -1
B -1 1
C 2 -1 -1
D 1
E 1
F 1 -1
Concentration changes
Stoichiometric matrix
Rate laws
dx N?v
dt N?v
Read the rate laws from the columns
Read the concentration changes from the rows
dA/dt -v1 - v2 dB/dt -v1 v4 dC/dt
2?v1 - v2 - v3 dD/dt v2 dE/dt v3 dF/dt
v3 - v4
vi(x,ei,ki)
x chemical species - concentrations v rate
laws k kinetic parameters N stoichiometric
matrix e catalysts (if any)
v v
v1 k1?A?B
v2 k2?A?C
v3 k3?C
v4 k4?F
E.g. dA/dt -k1?A?B - k2?A?C
65
From Chemical Reactions to Processes
r1 AB ?k1 CC r2 AC ?k2 D r3 C ?k3 EF r4 F
?k4 B
Write the coefficients by columns
Stoichiometric Matrix
interactions
N r1 r2 r3 r4
A -1 -1
B -1 1
C 2 -1 -1
D 1
E 1
F 1 -1
processes
For binary reactoins, first species in the column
does an input and produces result, second species
does an ouput, For unary reactions, species does
a tau action and produces result. No ternary
reactions.
A ?v1k1.(CC) ?v2k2.D ?a B !v1k1 ?b C
!v2k2 tk3(EF) ?c D 0 ?d E 0 ?e F
tk3.B ?f
Read the process interactions from the rows
Add a barb for counting and plotting
(Rate laws are implicit in stochastic semantics)
66
Stoichiometric Matrices Blow Up

• We can translate Chemistry to ODEs or Processes
• It is standard to go from chemical equations to
  ODEs via a stoichiometric matrix.
• It is similarly possible to go from chemical
  equations to processes via a stoichiometric
  matrix.
• But there is a better way
• Stoichiometric matrices blow-up exponentially for
  biochemical systems (unlike for ordinary chemical
  systems) because proteins have combinatorial
  state and complexed states are common.
• To avoid this explosion, we should describe
  biochemical systems compositionally without going
  through a stochiometric matrix (and hence without
  ODEs).

67
Complexes The ODE Way
A, B, C
ABC 69 ApBC ABC 69 ABpC ABC 69 ABCp ApBC 69
ApBpC ApBC 69 ApBCp ABpC 69 ApBpC ABpC 69
ABpCp ABCp 69 ApBCp ABCp 69 ABpCp ApBpC 69
ApBpCp ApBCp 69 ApBpCp ABpCp 69 ApBpCp
The matrix is very sparse, so the corresponding
ODE system is not dense. But it still has 2n
equations, one per species, plus conservation
equations (ABCApBCconstant, etc.).
domains
A 69 Ap B 69 Bp C 69 Cp
ABC ApBC ABpC ABCp ApBpC ApBCp ABpCp ApBpCp
2n
domainreactions
2n(2n-1)
ABC
System description is exponential in the number
of basic components.
reactions(twice number of edges in n-dim
hypercube)
species
complex
StoichiometricMatrix
N v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12 v13 v14 v15 v16 v17 v18 v19 v20 v21 v22 v23 v24
ABC
ApBC
ABpC
ABCp
ApBpC
ApBCp
ABpCp
ApBpCp
2n x 2n(2n-1)
(twice number of edges in n-dim hypercube, which
is given by Sloanes Sequence A001787
0,1,4,12,32,80,192,448,1024,2304,5120,11264,24576,
53248, 114688,245760,524288,1114112,2359296,498073
6,10485760, 22020096,46137344,96468992,201326592,4
19430400,872415232, 1811939328,3758096384,77846282
24 )
68
Complexes The Reactive System Way
A 69 Ap B 69 Bp C 69 Cp
A ?knAp Ap ?phA B ?knBp Bp ?phB C
?knCp Cp ?phC A B C
processes
domainreactions
System description is linear in the number of
basic components.
When the local domain reactions are not
independent, we can use lateral communication so
that each component is aware of the relevant
others.
(Its run-time behavior or analysis potentially
blows-up just as in the previous case, but its
description does not.)
69
Model Validation
70
Model Validation Simulation

• Basic stochastic algorithm Gillespie
• Exact (i.e. based on physics) stochastic
  simulation of chemical kinetics.
• Can compute concentrations and reaction times for
  biochemical networks.
• Stochastic Process Calculi
• BioSPi Shapiro, Regev, Priami, et. al.
• Stochastic process calculus based on Gillespie.
• BioAmbients Regev, Panina, Silverma, Cardelli,
  Shapiro
• Extension of BioSpi for membranes.
• Case study Lymphocytes in Inflamed Blood Vessels
  Lecaa, Priami, Quaglia
• Original analysis of lymphocyte rolling in blood
  vessels of different diameters.
• Case study Lambda Switch Celine Kuttler, IRI
  Lille
• Model of phage lambda genome (well-studied
  system).
• Case study VICE U. Pisa
• Minimal prokaryote genome (180 genes) and
  metabolism of whole VIrtual CEll, in stochastic
  p-calculus, simulated under stable conditions for
  40K transitions.
• Hybrid approaches
• Charon language UPenn
• Hybrid systems continuous differential equations
  discrete/stochastic mode switching.

71
Model Validation Program Analysis

• Causality Analysis
• Biochemical pathways, (concurrent traces such
  as the one here), are found in biology
  publications, summarizing known facts.
• This one, however, was automatically generated
  from a program written in BioSpi by comparing
  traces of all possible interactions. Curti,
  Priami, Degano, Baldari
• One can play with the program to investigate
  various hypotheses about the pathways.
• Control Flow Analysis
• Flow analysis techniques applied to process
  calculi.
• Overapproximation of behavior used to answer
  questions about what cannot happen.
• Analysis of positive feedback transcription
  regulation in BioAmbients Flemming Nielson.
• Probabilistic Abstract Interpretation
• DiPierro Wicklicky.

72
Model Validation Modelchecking

• Temporal
• Software verification of biomolecular systems (NA
  pump)Ciobanu
• Analysis of mammalian cell cycle (after Kohn) in
  CTL.Chabrier-Rivier Chiaverini Danos Fages
  Schachter
• E.g. is state S1 a necessary checkpoint for
  reaching state S2?
• Quantitative Simpathica/xssys Antioniotti Park
  Policriti Ugel Mishra
• Quantitative temporal logic queries of human
  Purine metabolism model.
• Stochastic Spring Parker Normal Kwiatkowska
• Designed for stochastic (computer) network
  analysis
• Discrete and Continuous Markov Processes.
• Process input language.
• Modelchecking of probabilistic queries.

Eventually(Always (PRPP 1.7 PRPP1)
implies steady_state() and
Eventually(Always(IMP lt 2 IMP1))
and Eventually(Always(hx_pool lt 10hx_pool1)))
73
What Reactive Systems Do For Us

• We can write things down precisely
• We can modularly describe high structural and
  combinatorial complexity (do programming).
• We can calculate and analyze
• Directly support simulation.
• Support analysis (e.g. control flow, causality,
  nondeterminism).
• Support state exploration (modelchecking).
• We can visualize
• Automata-like presentations.
• Petri-Net-like presentations.
• State Charts, Live Sequence Charts Harel
• Hierarchical automata.
• Scenario composition.

• We can reason
• Suitable equivalences on processes induce
  algebraic laws.
• We can relate different systems (e.g. equivalent
  behaviors).
• We can relate different abstraction levels.
• We can use equivalences for state minimization
  (symmetries).
• Disclaimers
• Some of these technologies are basically ready
  (medium-scale stochastic simulation and analysis,
  medium-scale nondeterministic and stochastic
  modelchecking).
• Others need to scale up significantly to be
  really useful. This is (has been) the challenge
  for computer scientists.

• Many approaches, same basic philosophy, tools
  being built
• ? Proc. Computational Methods in Systems Biology
  2003-2005

74
Conclusions
Q
The data are accumulating and the computers are
humming, what we are lacking are the words, the
grammar and the syntax of a new language D.
Bray (TIBS 22(9)325-326, 1997)
A

• The most advanced tools for computer process
  description seem to be also the best tools for
  the description of biomolecular systems.
  E.Shapiro (Lecture Notes)

75
References
MCB Molecular Cell Biology, Freeman. MBC
Molecular Biology of the Cell, Garland. Ptashne
A Genetic Switch. Davidson Genomic Regulatory
Systems. Milner Communicating and Mobile
Systems the Pi-Calculus. Regev Computational
Systems Biology A Calculus for Biomolecular
Knowledge (Ph.D. Thesis).
Papers BioAmbients a stochastic calculus with
compartments.Brane Calculi process calculi
with computation on the membranes, not inside
them. Bitonal Systems membrane reactions and
their connections to local patch
reactions. Abstract Machines of Systems
Biology the abstract machines implemented by
biochemical toolkits. www.luca.demon.co.uk/BioCom
puting.htm
76
(No Transcript)
77
(No Transcript)
78
Francois Hakim 3B
Parameters as in paper.Initially 1xA, 100xB, 3
copies of each gene.
Parameters as in paper.Initially 100xA, 1xB, 3
copies of each gene.Not really bistable.
79
François Hakim Fig415
80
François Hakim Fig415 in SPiM
( Francois and Hakim circuit of Fig 4 ) 5000.0
/ 2000 (--- ABCC ---) new ABCCltgt ( new
ABCCstayltgt ( !ABCC() ABCCstay()() (---
CCCC ---) new CCCCltgt ( new
tauDk_CCCC0.24ltgt ( new tf_cccc0.5ltltgtgt (
!CCCC() (new unb0.33ltgt (
tauDk_CCCC()() tf_ccccltunbgt unb()
CCCCltgt)) !tauDk_CCCCltgt (--- CC ---) new
CCltgt ( new tauDk_CC0.0014ltgt ( new
ABCCcpx0.74ltgt ( new CCCCcpx0.185ltgt ( 0.37/2
for homodimerization ) ( new CCBdk0.21ltgt (
!CC() tauDk_CC()() ABCCcpx()ABCCltgt
CCBdk()CCltgt CCCCcpx()CCCCltgt
CCCCcpxltgt() !tauDk_CCltgt (--- AB ---)
new ABltgt ( new tauDk_AB0.0057ltgt ( new
ABcpx0.066ltgt ( new ABCdk1.5ltgt ( !AB()
tauDk_AB()() ABCCcpxltgt() ABCdk()ABltgt
!tauDk_ABltgt (--- A ---) new Altgt ( new
tauDk_A0.00019ltgt ( !A() tauDk_A()()
ABcpx()ABltgt !tauDk_Altgt (--- B ---) new
Bltgt ( new tauDk_B0.077ltgt ( !B() tauDk_B()()
ABcpxltgt() CCBdkltgt() !tauDk_Bltgt (--- C
---) new Cltgt ( new tauDk_C0.023ltgt ( new
CCcpx0.07ltgt ( 0.14/2 for homodimerization
) ( !C() tauDk_C()() CCcpx()CCltgt
CCcpxltgt() ABCdkltgt() !tauDk_Cltgt
(--- a gate ---) new gateltltltgtgt,ltgt,ltgt,ltgtgt
( new gateBoundltltltgtgt,ltgt,ltgt,ltgt,ltgtgt (
!gate(a,b,tauCst,tauStim) a(a_)
(gateBoundlta,a_,b,tauCst,tauStimgt)
tauCst() (bltgt gatelta,b,tauCst,tauStimgt)
!gateBound(a,a_,b,tauCst,tauStim) a_ltgt
gatelta,b,tauCst,tauStimgt tauStim() (bltgt
gateBoundlta,a_,b,tauCst,tauStimgt) (--- gene a
--- ) new geneAltgt ( new geneAtauCst0.0000001
ltgt ( new geneAtauStim0.56ltgt ( !geneA()
gatelttf_cccc,A,geneAtauCst,geneAtauStimgt
!geneAtauCstltgt !geneAtauStimltgt (--- gene b
--- ) new geneBltgt ( new geneBtauCst0.43ltgt (
new geneBnoInput1.0ltltgtgt ( new
geneBnoStim1.0ltgt ( !geneB() geneBtauCst()
(Bltgt geneBltgt) !geneBtauCstltgt (--- gene c
--- ) new geneCltgt ( new geneCtauCst0.57ltgt (
new geneCnoInput1.0ltltgtgt ( new
geneCnoStim1.0ltgt ( !geneC() geneCtauCst()
(Cltgt geneCltgt) !geneCtauCstltgt geneAltgt
geneBltgt geneCltgt )))))))))) )))))))))) )))))))))
) ))))
81
(No Transcript)
82
Neg Gate Signal Response
neg(a,b) ?a th neg(a,b) te (tr(b)
neg(a,b))
inhibitory input
inhibition delay

back to initial state
constitutivetranscription
transcription factor
back to initial state
tr(b) !b tr(b) td
choice
binding
degradation
delay
h1 b100/a (at the fixpoint) matches Alons
numbersh0.01 b1/a is the self-feedback
instability pointih10 b900/a hence b
100h/a
a
83
Stochastic p-calculus

• Stochastic extension of p-calculus. C.Priami
• Associate a single parameter r (rate) in (0,
  infinity to each I/O activity a.
• The rate and the associated exponential
  distribution describes the stochastic behavior of
  the activity.
• a.P is replaced by a_at_r.P
• Exponential distribution
• guarantees the memoryless property the time at
  which a change of state occurs is independent of
  the time at which the last change of state
  occurred.
• Race condition
• is defined in a probabilistic competitive
  context all the activities that are enabled in a
  state compete and the fastest one
  (stochastically) succeeds.
• New implementation SPiM. A.Phillips. Paper at
  BioConcur