Luca Cardelli Microsoft Research Computability in Europe, Swansea 2006-07-05 www.luca.demon.co.uk

1
Luca CardelliMicrosoft ResearchComputabili
ty in Europe, Swansea2006-07-05www.luca.demon.c
o.uk
Artificial Biochemistry Biological Systems as
Reactive Systems
2
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
3
Stochastic Collectives
4
Stochastic Collectives

• Collective
• A large set of interacting finite state automata
• Not quite language automata (large set)
• Not quite cellular automata (interacting but
  not on a grid)
• Not quite process algebra (finite state and
  collective)
• Cf. multi-agent systems and swarm
  intelligence
• Stochastic
• Interactions have rates
• Not quite discrete (hundreds or thousands of
  components)
• Not quite continuous (non-trivial stochastic
  effects)
• Not quite hybrid (no switching between regimes)
• Very much like biochemistry
• Which is a large set of stochastically
  interacting molecules/proteins
• Are proteins finite state and subject to
  automata-like transitions?
• Lets say they are, at least because
• Much of the knowledge being accumulated in
  Systems Biology is described as state transition
  diagrams Kitano.

5
State Transitions
6
Even More State Transitions
http//www.expasy.ch/cgi-bin/show_thumbnails.pl
7
Interacting Automata
B1
A1
!a
_at_r1
?a
t_at_?1
t_at_?2
B2
B3
A2
t_at_?5
?b
_at_r2
!b
!c
_at_r3
A3
C1
C2
t_at_?3
?c
t_at_?4
C3
Communicating automata a graphical FSA-like
notation for finite state restriction-free
p-calculus processes. Interacting automata do
not even exchange values on communication. The
stochastic version has rates on communications,
and delays.
Finite state means no composition or
restriction inside recursion. Analyzable by
standard Markovian techniques, by first computing
the product automaton to obtain the underlying
finite Markov transition system. Buchholz
8
Interacting Automata Transition Rules
Exponential distribution P(delayt) 1-e -rt
Current State
Delay
Delay
Transition
t_at_r
t_at_r
r
a_at_r
a_at_r
Interaction
?a
!a
?a
!a
r
Q What kind of mass behavior can this
produce? (We need to understand that if want to
understand biochemical systems.)
9
Groupies and Celebrities
Groupie(wants to be like somebody different)
Celebrity(does not want to be like somebody else)
directive sample 0.1 200 directive plot A()
B() new a_at_1.0chan() new b_at_1.0chan() let A()
do !a A() or ?b B() and B() do !b B() or ?a
A() run 100 of (A() B())
directive sample 0.1 200 directive plot A()
B() new a_at_1.0chan() new b_at_1.0chan() let A()
do !a A() or ?a B() and B() do !b B() or ?b
A() run 100 of (A() B())
A stochastic collective of celebrities
A stochastic collective of groupies
always eventually deadlock
always equilibrium
Unstable because within an A majority, an A has
difficulty finding a B to emulate, but the few
Bs have plenty of As to emulate, so the
majority may switch to B. Leads to deadlock when
everybody is in the same state and there is
nobody different to emulate.
Stable because as soon as a A finds itself in the
majority, it is more likely to find somebody in
the same state, and hence change, so the majority
is weakened.
10
Both Together
A way to break the deadlocks Groupies with just
a few Celebrities
!a
!a
directive sample 10.0 1000 directive plot Ga()
Gb() Ca() Cb() new a_at_1.0chan() new
b_at_1.0chan() let Ca() do !a Ca() or ?a
Cb() and Cb() do !b Cb() or ?b Ca() let Ga()
do !a Ga() or ?b Gb() and Gb() do !b Gb()
or ?a Ga() run 1 of (Ca() Cb()) run 100
of (Ga() Gb())
Ca
Ga
A few (2)Celebrities
ManyGroupies
?b
?a
?a
?b
Cb
Gb
!b
!b
never deadlock
11
Hysteric Groupies
Regularity can arise not far from chaos
We can get more regular behavior from groupies if
they need more convincing, or hysteresis
(history-dependence), to switch states.
!a
directive sample 10.0 1000 directive plot Ga()
Gb() new a_at_1.0chan() new b_at_1.0chan() let Ga()
do !a Ga() or ?b ?b Gb() and Gb() do !b
Gb() or ?a ?a Ga() let Da() !a Da() and
Db() !b Db() run 100 of (Ga() Gb()) run 1
of (Da() Db())
a solid threshold to observe switching
A
?a
?b
?a
?b
B
!b
!a
!b
(With doping to break deadlocks)
N.B. It will not oscillate without doping (noise)
regular oscillation
!a
directive sample 10.0 1000 directive plot Ga()
Gb() new a_at_1.0chan() new b_at_1.0chan() let Ga()
do !a Ga() or ?b ?b ?b Gb() and Gb() do
!b Gb() or ?a ?a ?a Ga() let Da() !a
Da() and Db() !b Db() run 100 of (Ga()
Gb()) run 1 of (Da() Db())
A
?a
?b
?a
?b
?a
?b
B
!b
12
Hysteric 3-Way Groupies
directive sample 3.0 1000 directive plot A()
B() C() new a_at_1.0chan() new b_at_1.0chan() new
c_at_1.0chan() let A() do !a A() or ?c ?c
C() and B() do !b B() or ?a ?a A() and C()
do !c C() or ?b ?b B() let Da() !a
Da() and Db() !b Db() and Dc() !c Dc() run
100 of (A() B() C()) run 1 of (Da() Db()
Dc())
!a
!b
?a
?a
A
B
!a
!b
?c
?b
?c
?b
Da
Db
Dc
C
(Still with doping)
!c
!c
N.B. It will not oscillate without doping (noise)
13
Oscillation as Emergence
Just 2 of the hysteric groupies do not oscillate
regularly at all!
Nor 16
Dotted lines indicate cross sections where one
may look for evidence of alternation.
new a_at_1.0chan() new b_at_1.0chan() let A() do
!a A() or ?b ?b ?b B() and B() do !b B()
or ?a ?a ?a A() let As() !a As() and Bs()
!b Bs() run 64 of (A() B()) run 1 of
(As() Bs())
Pretty good with 64
14
Distributions can be Programmed
Exercise (hard) Build a small automaton where
one state has an occupation distribution like
this
Or, more specifically, build a 3-state, A-B-C,
automaton such that
B B(A-C)
15
Semantics of Collective Behavior
16
Micromodels Continuous Time Markov Chains

• The underlying semantics of stochastic p-calculus
  (and stochastic interacting automata). Well
  established in many ways.
• Automata with rates on transitions.
• The correct semantics for chemistry,
  executable.
• Gillespie stochastic simulation algorithm
• Lots of advantages
• Compositional, compact, mechanistic, etc.
• But do not give a good sense of collective
  properties.
• Yes one can do simulation.
• Yes one can do program analysis.
• Yes one can do modelchecking.
• But somewhat lacking in analytical properties
  and predictive power.

17
Macromodels Ordinary Differential Equations

• They always ask
• Yes, but how does you automata model relate to
  the 75 ODE models in the literature?
• Going from processes/automata to ODEs directly
• In principle just write down the Rate Equation
  Calder, Hillston
• - Determine the set of all possible states S of
  each process.
• - Determine the rates of the transitions between
  such states.
• - Let S be the number of processes in state S
  as a function of time.
• - Define for each state S
• S (rate of change of the number of
  processes in state S)
• Cumulative rate of transitions from any state
  S to state S, times S,
• minus cumulative rate of transitions from S to
  any state S, times S.
• - Intuitive (rate inflow minus outflow), but
  often clumsy to write down precisely.
• But why go to the trouble?
• If we first convert processes to chemical
  reactions, then we can convert to ODEs by
  standard means!

18
Macromodel of InteractionLaw of Mass Interaction
directive sample 0.01 1000 directive plot D()
A1() A2() A4() A8() new c1_at_1.0 chan() new
c2_at_2.0 chan() new c4_at_4.0 chan() new c8_at_8.0
chan() let D() delay_at_1.0 let A1() ?c1 and
B1() !c1 let A2() ?c2 and B2() !c2 let A4()
?c4 and B4() !c4 let A8() ?c8 and B8() !c8
?1,2,4,8 A0100
Slow
Fast
?1
A01000
decay
N interactions on the same channel are faster
than N delays at the same rate (on N independent
channels).
interaction
?1,2,4,8
Each C has ? chances to decay per second (no
matter how many other Cs there are), but each A
has B? chances to interact per second it
depends on how many Bs there are.
directive sample 0.01 1000 directive plot D()
A1() A2() A4() A8() new c1_at_1.0 chan() new
c2_at_2.0 chan() new c4_at_4.0 chan() new c8_at_8.0
chan() let D() delay_at_1.0 let A1() ?c1 and
B1() !c1 let A2() ?c2 and B2() !c2 let A4()
?c4 and B4() !c4 let A8() ?c8 and B8()
!c8 run 1000 of (D() A1() B1() A2()
B2() A4() B4() A8() B8())
19
From Chemistry to ODEs
20
Chemical Reactions
A r B1 Bn A1 A2 r B1 Bn A
A r B1 Bn
Degradation A -rA Asymmetric
Collision Ai -rA1A2 Symmetric
Collision A -rA(A-1)
Exponential Decay Mass Action Law Mass Action
Law
(assuming A?Bi?Aj for all i,j)
No other reactions!
Chapter IV Chemical Kinetics David A.
Reckhow , CEE 572 Course ... reactions may be
either elementary or non-elementary. Elementary
reactions are those reactions that occur exactly
as they are written, without any intermediate
steps. These reactions almost always involve just
one or two reactants. ... Non-elementary
reactions involve a series of two or more
elementary reactions. Many complex environmental
reactions are non-elementary. In general,
reactions with an overall reaction order greater
than two, or reactions with some non-integer
reaction order are non-elementary.
THE COLLISION THEORY OF REACTION RATES
www.chemguide.co.uk The chances of all this
happening if your reaction needed a collision
involving more than 2 particles are remote. All
three (or more) particles would have to arrive at
exactly the same point in space at the same time,
with everything lined up exactly right, and
having enough energy to react. That's not likely
to happen very often!
21
From Reactions to ODEs
CAVEAT A deterministic approximation of a
stochastic system (i.e. possibly misleading)
v1 AB ?k1 CC v2 AC ?k2 D v3 C ?k3 EF v4
FF ?k4 B
Write the coefficients by columns
Stoichiometric Matrix
N v1 v2 v3 v4
A -1 -1
B -1 1
C 2 -1 -1
D 1
E 1
F 1 -2
Quantity changes
Stoichiometric matrix
Rate laws
X N?l
Set a rate law for each reaction
(Degradation/Asymmetric/Symmetric)
Read the concentration changes from the rows
A -l1 - l2 B -l1 l4 C 2l1 - l2 -
l3 D l2 E l3 F l3 - 2l4
X chemical species - quantity of molecules l
rate laws k kinetic parameters N stoichiometric
matrix
l l
l1 k1AB
l2 k2AC
l3 k3C
l4 k4F(F-1)/2
E.g. A -k1AB - k2AC
22
From Processes to Chemistry
23
Chemical Ground Form (CGF)
E X1M1, , XnMn Definitions (n0) M
p1P1 Å Å pnPn Molecules (n0) P
X1 Xn Solutions (n0) p
tr ?n(r) !n(r) Interactions (delay, input,
output) CGF E,P Definitions with
Initial Conditions
Å is stochastic choice (vs. for chemical
reactions) 0 is the null solution (P0 0P P)
and null molecule (MÅ0 0ÅM M) (t0P 0)
Xi are distinct in E Each name n is assigned a
fixed rate r n(r)
(To translate chemistry back to processes we need
a bit more than simple automata we may have
on the right of , that is we may need after
p.)
Ex interacting automata (which are CGFs using
only in initial conditions)
Automaton in state A
A !aA Å ?bB B !bB Å ?aA AABB
Automaton in state B
Initial conditions 2A and 2B
24
CGF Semantics
Reduction
E, (X1 P) r E, (P1 P) if E º
X1trP1ÅM1, E E, (X1 X2 P) r E, (P1 P2
P) if E º X1?n(r)P1ÅM1, E1 º
X2!n(r)P2ÅM2, E2 E,P r E,P if E,P º
E,P1 Ù E,P1 r E,P2 Ù E,P2 º E,P Þ
Structural Congruence
º is an equivalence relation E, E º E, E M
Å M º M Å M P P º P P
E º E Ù P º P Þ E,P º E,P E º E Ù M º M
Þ XM, E º XM, E M º M Ù P º P Þ pP Å
M º pP Å M P º P Þ X P º X P
E (A !aA Å ?bB B !bB Å ?aA) E,
(ABB) r(a) E, (AAB) r(b) E, (ABB)
r(b) E, (BBB)
25
Automata to Chemistry
AB BB BA AA
ABd BBd BAd AAd
AC CC CB BB BA AA
ACd CCd CBd BBd BAd AAd
26
Three Main Cases
Unary reactions. These are not finite state
systems, but finite species systems are ok!
E X tr(X X)
C(E) X r X X
Unbounded state, but only 1 species.No problem!
Binary reactions. The same interaction can occur
multiple times and must be taken into account
E A ?nB Å ?nB C !nD
C(E) A C r(n) B D A C r(n) B D
That is A C 2r(n) B D
Symmetric reactions
E X !a0 Å ?aY
C(E) X X 2r(a) Y The rate of a was
pre-halved and must be restored.
27
CGF to Chemistry
E X1M1, , XnMn Definitions (n0) M
p1P1 Å Å pnPn Molecules (n0) P
X1 Xn Solutions (n0) p
tr ?n(r) !n(r) Interactions (delay, input,
output) CGF E,P Definitions with
Initial Conditions
Each X in E is seen as a separate species.
Chemical reactions for E (N.B. m is a
multiset, and P is P with all the changed to )
ChG(E) (X r P) s.t. (X º trP Å ) Î Em
Èm (X Y r P Q) s.t. X?Y, á(X º ?n(r)P Å
),(Y º !n(r)Q Å )ñ Î E2m Èm (X X 2r P
Q) s.t. (X º ?n(r)P Å º !n(r)Q Å )ñ Î Em
Initial conditions for P
ChG(P) P
28
From Processes to ODEs
29
Nonlinear Transitions
30
Basic Nonlinear Transition
!c
_at_s
?c
directive sample 0.02 1000 directive plot B()
A() val s1.0 new c_at_schan let A() ?c
B() and B() !cB() run (1000 of A() 1 of B())
A
B
N.B. needs at least 1 B to get started.
A ?c(s)B B !c(s)B
SPiM
AB s BB
A -sABB sAB
Matlab continuous_sys_generator
interval/step 00.0010.0 (A) dx1/dt -
x1x2 1000.0 (B) dx2/dt x1x2 1.0
31
Bell Exercise
Build a small network where one node has a
distribution like B()
A ?b(1)B B !b(1)B Å ?c(1)C C !c(1)C
SPiM
B B(A-C)
AB 1 BB BC 1 CC
directive sample 0.0025 1000 directive plot B()
A() C() new b_at_1.0chan new c_at_1.0chan let A()
?b B() and B() do !bB() or ?c C() and C()
!cC() run ((10000 of A()) B() C())
A -ABB AB-BC C BC
Matlab continuous_sys_generator
interval/step 00.0000010.0025 (A) dx1/dt
-x1x2 10000.0 (B) dx2/dt x1x2
x2x3 1.0 (C) dx3/dt x2x3 1.0
32
Oscillator
directive sample 0.1 1000 directive plot A1()
A2() A3() val r1.0 val s1.0 new a1_at_schan
new a2_at_schan new a3_at_schan let A1() do
!a1A1() or delay_at_rA2() or ?a2 ?a2 A2() and
A2() do !a2A2() or delay_at_rA3() or ?a3 ?a3
A3() and A3() do !a3A3() or delay_at_rA1() or
?a1 ?a1 A1() run 1000 of A1()
Sustained Determinisitic Oscillation
N.B. this does not deadlock!
Robust Stochastic Oscillation
Matlab continuous_sys_generator
SPiM
Matlab continuous_sys_generator
interval/step 00.0012.0 (A) dx1/dt -x1 -
x1x2 x3 x6x1 100.0 (B) dx2/dt -x2 - x2x3
x1 x4x2 0.0 (C) dx3/dt -x3 - x3x1 x2
x5x3 0.0 (A) dx4/dt -x4x2
x1x2 0.0 (B) dx5/dt -x5x3 x2x3
0.0 (C) dx6/dt -x6x1 x3x1 0.0
SPiM
33
Epidemics
Kermack, W. O. and McKendrick, A. G. "A
Contribution to the Mathematical Theory of
Epidemics." Proc. Roy. Soc. Lond. A 115, 700-721,
1927. http//mathworld.wolfram.com/Kermack-McKend
rickModel.html
34
Epidemics
directive sample 500.0 1000 directive plot
Recovered() Susceptible() Infected() new
infect _at_0.001chan() val recover 0.03 let
Recovered() ?infect Recovered() and
Susceptible() ?infect Infected() and
Infected() do !infect Infected() or
?infect Infected() or delay_at_recover
Recovered() run (200 of Susceptible() 2 of
Infected())
35
ODE
lt?xml version"1.0" encoding"UTF-8"?gt lt!--
Created by SBML API 2.0(a17.0) --gt ltsbml
level"2" version"1" xmlns"http//www.sbml.org/s
bml/level2" xmlnscelldesigner"http//www.sbml.or
g/2001/ns/celldesigner"gt ltmodel
id"test"gt ltannotationgt ltcelldesignermodelVersion
gt2.5lt/celldesignermodelVersiongt ltcelldesignermod
elDisplay sizeX"600" sizeY"400"/gt ltcelldesigner
listOfCompartmentAliases/gt ltcelldesignerlistOfCom
plexSpeciesAliases/gt ltcelldesignerlistOfSpeciesAl
iasesgt ltcelldesignerspeciesAlias id"sa18"
species"s9"gt ltcelldesigneractivitygtinactivelt/cel
ldesigneractivitygt ltcelldesignerbounds h"25.0"
w"70.0" x"36.0" y"152.5"/gt ltcelldesignerview
state"usual"/gt ltcelldesignerusualViewgt ltcelldesi
gnerinnerPosition x"0.0" y"0.0"/gt ltcelldesigner
boxSize height"25.0" width"70.0"/gt ltcelldesigne
rsingleLine width"1.0"/gt ltcelldesignerpaint
color"ffccff66" scheme"Color"/gt lt/celldesigneru
sualViewgt ltcelldesignerbriefViewgt ltcelldesigneri
nnerPosition x"0.0" y"0.0"/gt ltcelldesignerboxSi
ze height"60.0" width"80.0"/gt ltcelldesignersing
leLine width"0.0"/gt ltcelldesignerpaint
color"3fff0000" scheme"Color"/gt lt/celldesignerb
riefViewgt lt/celldesignerspeciesAliasgt ltcelldesign
erspeciesAlias id"sa19" species"s10"gt ltcelldesi
gneractivitygtinactivelt/celldesigneractivitygt ltce
lldesignerbounds h"25.0" w"70.0" x"36.0"
y"58.5"/gt ltcelldesignerview state"usual"/gt ltcel
ldesignerusualViewgt ltcelldesignerinnerPosition
x"0.0" y"0.0"/gt ltcelldesignerboxSize
height"25.0" width"70.0"/gt ltcelldesignersingleL
ine width"1.0"/gt ltcelldesignerpaint
color"ffccff66" scheme"Color"/gt lt/celldesigneru
sualViewgt ltcelldesignerbriefViewgt ltcelldesigneri
nnerPosition x"0.0" y"0.0"/gt ltcelldesignerboxSi
ze height"60.0" width"80.0"/gt ltcelldesignersing
leLine width"0.0"/gt ltcelldesignerpaint
color"3fff0000" scheme"Color"/gt lt/celldesignerb
riefViewgt lt/celldesignerspeciesAliasgt ltcelldesign
erspeciesAlias id"sa20" species"s11"gt ltcelldesi
gneractivitygtinactivelt/celldesigneractivitygt ltce
lldesignerbounds h"25.0" w"70.0" x"273.0"
y"40.5"/gt ltcelldesignerview state"usual"/gt ltcel
ldesignerusualViewgt ltcelldesignerinnerPosition
x"0.0" y"0.0"/gt ltcelldesignerboxSize
height"25.0" width"70.0"/gt ltcelldesignersingleL
ine width"1.0"/gt ltcelldesignerpaint
color"ffccff66" scheme"Color"/gt lt/celldesigneru
sualViewgt ltcelldesignerbriefViewgt ltcelldesigneri
nnerPosition x"0.0" y"0.0"/gt ltcelldesignerboxSi
ze height"60.0" width"80.0"/gt ltcelldesignersing
leLine width"0.0"/gt ltcelldesignerpaint
color"3fff0000" scheme"Color"/gt lt/celldesignerb
riefViewgt lt/celldesignerspeciesAliasgt lt/celldesig
nerlistOfSpeciesAliasesgt ltcelldesignerlistOfGrou
ps/gt ltcelldesignerlistOfProteins/gt ltcelldesigner
listOfGenes/gt ltcelldesignerlistOfRNAs/gt ltcelldesi
gnerlistOfAntisenseRNAs/gt ltcelldesignerlistOfBlo
ckDiagrams/gt lt/annotationgt ltlistOfCompartmentsgt ltc
ompartment id"default"/gt lt/listOfCompartmentsgt ltl
istOfSpeciesgt ltspecies compartment"default"
id"s9" initialAmount"200.0" name"S"gt ltannotatio
ngt ltcelldesignerpositionToCompartmentgtinsidelt/cel
ldesignerpositionToCompartmentgt ltcelldesignerspe
ciesIdentitygt ltcelldesignerclassgtSIMPLE_MOLECULElt
/celldesignerclassgt ltcelldesignernamegtSlt/celldes
ignernamegt lt/celldesignerspeciesIdentitygt lt/anno
tationgt lt/speciesgt ltspecies compartment"default"
id"s10" initialAmount"2.0" name"I"gt ltannotation
gt ltcelldesignerpositionToCompartmentgtinsidelt/cell
designerpositionToCompartmentgt ltcelldesignerspec
iesIdentitygt ltcelldesignerclassgtSIMPLE_MOLECULElt/
celldesignerclassgt ltcelldesignernamegtIlt/celldesi
gnernamegt lt/celldesignerspeciesIdentitygt ltcellde
signerlistOfCatalyzedReactionsgt ltcelldesignercat
alyzed reaction"re26"/gt lt/celldesignerlistOfCata
lyzedReactionsgt lt/annotationgt lt/speciesgt ltspecies
compartment"default" id"s11" initialAmount"0.0"
name"R"gt ltannotationgt ltcelldesignerpositionToCo
mpartmentgtinsidelt/celldesignerpositionToCompartme
ntgt ltcelldesignerspeciesIdentitygt ltcelldesignerc
lassgtSIMPLE_MOLECULElt/celldesignerclassgt ltcelldes
ignernamegtRlt/celldesignernamegt lt/celldesignersp
eciesIdentitygt lt/annotationgt lt/speciesgt lt/listOfSp
eciesgt ltlistOfReactionsgt ltreaction id"re25"
reversible"false"gt ltannotationgt ltcelldesignerrea
ctionTypegtSTATE_TRANSITIONlt/celldesignerreactionT
ypegt ltcelldesignerbaseReactantsgt ltcelldesignerba
seReactant alias"sa19" species"s10"gt ltcelldesign
erlinkAnchor position"N"/gt lt/celldesignerbaseRe
actantgt lt/celldesignerbaseReactantsgt ltcelldesigne
rbaseProductsgt ltcelldesignerbaseProduct
alias"sa20" species"s11"gt ltcelldesignerlinkAnch
or position"S"/gt lt/celldesignerbaseProductgt lt/ce
lldesignerbaseProductsgt ltcelldesignerconnectSche
me connectPolicy"direct"gt ltcelldesignerlistOfLin
eDirectiongt ltcelldesignerlineDirection index"0"
value"unknown"/gt lt/celldesignerlistOfLineDirecti
ongt lt/celldesignerconnectSchemegt ltcelldesignerli
ne color"ff000000" width"1.0"/gt lt/annotationgt ltl
istOfReactantsgt ltspeciesReference
species"s10"gt ltannotationgt ltcelldesigneraliasgtsa
19lt/celldesigneraliasgt lt/annotationgt lt/speciesRef
erencegt lt/listOfReactantsgt ltlistOfProductsgt ltspeci
esReference species"s11"gt ltannotationgt ltcelldesig
neraliasgtsa20lt/celldesigneraliasgt lt/annotationgt
lt/speciesReferencegt lt/listOfProductsgt ltkineticLawgt
ltmath xmlns"http//www.w3.org/1998/Math/MathML"gt
ltapplygt lttimes/gt ltcngt0.03lt/cngt ltcigts10lt/cigt lt/app
lygt lt/mathgt lt/kineticLawgt lt/reactiongt ltreaction
id"re26" reversible"false"gt ltannotationgt ltcellde
signerreactionTypegtSTATE_TRANSITIONlt/celldesigner
reactionTypegt ltcelldesignerbaseReactantsgt ltcelld
esignerbaseReactant alias"sa18"
species"s9"gt ltcelldesignerlinkAnchor
position"N"/gt lt/celldesignerbaseReactantgt lt/cell
designerbaseReactantsgt ltcelldesignerbaseProducts
gt ltcelldesignerbaseProduct alias"sa19"
species"s10"gt ltcelldesignerlinkAnchor
position"S"/gt lt/celldesignerbaseProductgt lt/celld
esignerbaseProductsgt ltcelldesignerconnectScheme
connectPolicy"direct"gt ltcelldesignerlistOfLineDi
rectiongt ltcelldesignerlineDirection index"0"
value"unknown"/gt lt/celldesignerlistOfLineDirecti
ongt lt/celldesignerconnectSchemegt ltcelldesignerli
ne color"ff000000" width"1.0"/gt ltcelldesignerli
stOfModificationgt ltcelldesignermodification
aliases"sa19" modifiers"s10" targetLineIndex"-1
,0" type"STATE_TRANSITION"gt ltcelldesignerconnect
Scheme connectPolicy"direct"gt ltcelldesignerlistO
fLineDirectiongt ltcelldesignerlineDirection
index"0" value"unknown"/gt lt/celldesignerlistOfL
ineDirectiongt lt/celldesignerconnectSchemegt ltcelld
esignerlinkTarget alias"sa19"
species"s10"gt ltcelldesignerlinkAnchor
position"SW"/gt lt/celldesignerlinkTargetgt ltcellde
signerline color"ff000000" width"1.0"/gt lt/celld
esignermodificationgt lt/celldesignerlistOfModific
ationgt lt/annotationgt ltlistOfReactantsgt ltspeciesRef
erence species"s9"gt ltannotationgt ltcelldesigneral
iasgtsa18lt/celldesigneraliasgt lt/annotationgt lt/spec
iesReferencegt lt/listOfReactantsgt ltlistOfProductsgt
ltspeciesReference species"s10"gt ltannotationgt ltcel
ldesigneraliasgtsa19lt/celldesigneraliasgt lt/annota
tiongt lt/speciesReferencegt lt/listOfProductsgt ltlistO
fModifiersgt ltmodifierSpeciesReference
species"s10"gt ltannotationgt ltcelldesigneraliasgtsa
19lt/celldesigneraliasgt lt/annotationgt lt/modifierSp
eciesReferencegt lt/listOfModifiersgt ltkineticLawgt ltm
ath xmlns"http//www.w3.org/1998/Math/MathML"gt lta
pplygt lttimes/gt ltcngt0.001lt/cngt ltcigts9lt/cigt ltcigts10lt
/cigt lt/applygt lt/mathgt lt/kineticLawgt lt/reactiongt lt/
listOfReactionsgt lt/modelgt lt/sbmlgt
S ?i(t)I I !i(t)I Å ?i(t)I Å trR R
?i(t)R

• Cell Designer ODE Solver output for reactions
• S I t I I
• I r R
• with t 0.001 r 0.03 S200 I2

S I t I I I I t I I I r R R I t R
I
useless
useless reactions
S -tSII tSI-rI R rI

• ODE Solver output for
• S -tSI
• I tSI-rI
• R rI
• with t 0.001 r 0.03 S0200 I02

Automata match the standard ODE model!
Matlab continuous_sys_generator
36
Simplified Model
not useless!
S ?i(t)I I !i(t)I Å trR R 0
?infect
!infect
Susceptible
useless
useless
Infected
?infect
_at_recover
Not totally obvious that one could have
simplified the automata model.
S I t I I I r R
Recovered
?infect
S -tSII tSI-rI R rI
directive sample 500.0 1000 directive plot
Recovered() Susceptible() Infected() new
infect _at_0.001chan() val recover 0.03 let
Recovered() () and Susceptible()
?infect Infected() and Infected() do
!infect Infected() or delay_at_recover
Recovered() run (200 of Susceptible() 2 of
Infected())
Same ODE, hence equivalent automata models.
37
Lotka-Volterra
38
Predator-Prey
directive sample 5000.0 1000 directive plot
Carnivor() Herbivor() val mortality 0.01 val
breeding 0.01 val predation 0.01 new cull
_at_predationchan() let Herbivor() do
delay_at_breeding (Herbivor() Herbivor()) or
?cull () and Carnivor() do
delay_at_mortality () or !cull (Carnivor()
Carnivor()) run replicate delay_at_0.01
(Herbivor() Carnivor())
_at_breeding
Herbivor
?cull
_at_predation
!cull
Carnivor
_at_mortality
Since predator and prey drive each other to
exinction (stochastically), we restart the
popolations periodically.
An unbounded state system!
(This is a case where the continuous system
oscillates and the stochastic one does not! We
have seen examples of the opposite situation.)
39
ODE
H tb (HH) Å ?c(p)0 C tm0 Å !c(p)(CC)
H b H H C m 0 H C p C C
H bH-pHCC -mCpHC
Automata match the Lotka-Volterra model (with BD)
40
Laws by ODEs
41
Choice Law by ODEs
t?B Å tµB t?µB
t?
t?µ
A
A
tµ
A t?µB
A t?B Å tµB
A ?µ B
A ? B A µ B
A -(?µ)A B (?µ)A
A -?A - µA B ?A µA

42
Idle Delay Law by ODEs
A t?A Å tµB A tµB
t?
A
B
A
B
tµ
tµ
A t?A Å tµB
A tµB
A µ B
A ? A A µ B
A -µAB µA
A -µAB µA

43
Idle Interaction Law by ODEs
!c
?c
!c
?c
?c
A
B
C
A
B
C
It may seem like A should decrease half as fast,
but NO! Two ways to explain -State A is
memoryless of any past idling. - Activity on
c is double
A ?cB C !cC
A ?cA Å ?cB C !cC
directive sample 6.0 1000 directive plot A() new
c_at_1.0chan let A() ?c B() and B() () and
C() !c C() run (C() 1000 of A())
directive sample 6.0 1000 directive plot A() new
c_at_1.0chan let A() do ?c B() or ?c A() and
B() () and C() !c C() run (C() 1000 of
A())
AC r BC
AC r AC AC r BC
A -rACB rAC C 0
A -rACB rAC C 0
44
Asynchronous Interleaving
Hermanns Interactive Markov Chains. Sec 4.1.2
t?B tµD t?(B tµD) tµ(t?B D)
A10
directive sample 4.0 10000 directive plot A()
B() C() D() let A() delay_at_1.0 B() and B()
() let C() delay_at_2.0 D() and D() () run
1000 of (A() C())
_at_1
1000
1000
A
B
C10
_at_2
1000
1000
C
D
directive sample 4.0 10000 directive plot
?YA B() ?YC D() Y() A() C() new YA_at_1.0chan
new YC_at_1.0chan let A() do delay_at_1.0 B() or
?YA and B() () let C() do delay_at_2.0 D() or
?YC and D() () let Y() do delay_at_1.0
(B() C()) or delay_at_2.0 (A() D()) or
?YA or ?YC run 1000 of Y()
333
_at_2
333
_at_1
Y0
D
1000
Y
333
667
667
YA
B
_at_2
YC
_at_1
667
Amazingly, the Bs and the Ds from the two
branches sum up to exponential distributions
45
Asynchronous Interleaving Law by ODEs
t?B tµD t?(B tµD) tµ(t?B D)
Want to show that B and D on both sides have the
same behavior (equal quantities of B and D
produced at all times)
A1 t?B C1 tµD A1 C1
Y t?(B C2) Å tµ(A2 D) C2 tµD A2
t?B Y
B C2 - µC2 ?A2 C2 - D
µE2 ?A2
Y ? B C2 Y µ A2 D C2 µ D A2 ? B Y
A1 ? B C1 µ D A1 C1
Y -?Y-µY A2 µY-?A2 B
?Y?A2 C2 ?Y-µC2 D µYµC2
A1 -?A1 B ?A1 C1 -µC1 D
µC1
YA2 YA2 -?Y-µYµY-?A2
-?Y-?A2 -?YA2 YA2 decays
exponentially!
YA2 -?YA2 B ?YA2 YC2
-µYC2 D µYC2
E2 (C2 µC2)/? E2 (A2
?A2)/µ B ?(A2 ?A2)/µ ?A2
?
B and D have equal time evolutions on the two
sides provided that A1YA2 and C1YC2.
This imposes the constraint, in particular, that
A10YA20 and C10YC20 (at time
zero). The initial conditions of the right hand
system specify that A20C200 (since only Y
is present). Therefore, we obtain that
A10C10Y0.
B E1 - A1 -E1 -
A1 D -E1 - C1
B ?E2 µE2 - A2 (?µ)E2 -
A2 - E2- A2 D - E2-
C2
Equal ODEs up to a change of variables (A1YA2,
C1YC2), where B and D have the same
behavior! They are related in the same way.
So, for example, if we run a stochastic
simulation of the left hand side with 1000A1 and
1000C1, we obtain the same curves for B and D
than a stochastic simulation of the right hand
side with 1000Y.
46
Parametric Form
47
Chemical Parametric Form (CPF)
E X1(p1)M1, , Xn(pn)Mn
Definitions (n 0) M p1P1 Å Å pnPn
Molecules (n 0) P
X1(p1) Xn(pn) Solutions (n
0) p tr ?n(p) !n(p)
Interactions CPF E,P with
initial conditions
Å is stochastic choice (vs. for chemical
reactions) 0 is the null solution (P0 0P P)
and null molecule (MÅ0 0ÅM M) (t0P 0)
Xi are distinct in E, p are vectors of names p
are vectors of distinct names when in binding
position Each free name n in E is assigned a
fixed rate r written either n(r), or rCPF(n)r.
A translation from CPF to CGF exists (expanding
all possible instantiation of parameters from the
initial conditions) An incremental translation
algorithm exists (expanding on demand from
initial conditions)
48
Repressilator ODEs
simplifying (N is the quantity of each of the 3
gates)
Neg(a,b) ?a Inh(a,b) Å te (Tr(b)
Neg(a,b)) Inh(a,b) th Neg(a,b) Tr(b) !b
Tr(b) Å tg 0 Neg(x(r),y(r)) Neg(y(r),z(r))
Neg(z(r),x(r))
Neg/x,y hN (hrTr/x)Neg/x,y Neg/y,z
hN (hrTr/y)Neg/y,zNeg/z,x hN
(hrTr/z)Neg/z,xTr/x eNeg/z,x -
gTr/x Tr/y eNeg/x,y - gTr/y Tr/z
eNeg/y,z - gTr/z
Neg/x,y e Tr/y Neg/x,y Neg/y,z e Tr/z
Neg/y,z Neg/z,x e Tr/x Neg/z,x Tr/x
Neg/x,y r Tr/x Inh/x,y Tr/y Neg/y,z r Tr/y
Inh/y,z Tr/z Neg/z,x r Tr/z Inh/z,x
Inh/x,y h Neg/x,y Inh/y,z h Neg/y,z Inh/z,x h
Neg/z,x Tr/x g 0 Tr/y g 0 Tr/z g 0 Neg/x,y
Neg/y,z Neg/z,x
No sustained oscillations (with SPiM
parameters). But see ElowitzLeibler.
Neg/x,y -rTr/xNeg/x,y
hInh/x,yNeg/y,z -rTr/yNeg/y,z
hInh/y,zNeg/z,x -rTr/zNeg/z,x
hInh/z,xInh/x,y rTr/xNeg/x,y -
hInh/x,y Inh/y,z rTr/yNeg/y,z -
hInh/y,z Inh/z,x rTr/zNeg/z,x -
hInh/z,x Tr/x eNeg/z,x - gTr/x Tr/y
eNeg/x,y - gTr/y Tr/z eNeg/y,z -
gTr/z
49
Groupies ODE
50
Doped Groupies ODE
Q What does this do?
A !a(r)A Å ?b(r)B B !b(r)B Å ?a(r)A
Ad !a(r)Ad Bd !b(r)Bd
directive sample 10.0 1000 directive plot Ga()
Gb() Da() Db() new a_at_1.0chan() new
b_at_1.0chan() let Ga() do !a Ga() or ?b
Gb() and Gb() do !b Gb() or ?a Ga() let Da()
!a Da() and Db() !b Db() run 1 of (Da()
Db()) run 100 of (Ga() Gb())
AB r AA BA r BB
ABd r BBd BAd r AAd
A rAB-rBA-rABdrBAd B
rBA-rAB-rBAdrABd
Ad 0 Bd 0
Ad,Bd are constant assume them both k
At B0 A-rkA, BrkA At
AB AB0 At AB AB0
A -rk(A-B) B rk(A-B)
Stochastic Answer bounded random walk
Deterministic Answer convergence and stability
-1.0(y1-y2), 1.0 1.0(y1-y2), 0.0
00.0033.0 r1.0 k1.0 dx1/dt -(x1-x2),
200.0 dx2/dt (x1-x2), 0.0
ODE predicts converging stable equilibrium at
AB instead of the total chaos observed in
the stochastic system!
For k0 (no dope), predicts deadlock AB0
but at any value of A, which is definitely not
true in the stochastic system.
Matlab continuous_sys_generator
DESSolver
SPiM
51
Conclusions
52
Conclusions

• Stochastic Collectives
• Complex global behavior from simple components
• Emergence of collective functionality from
  non-functional components
• (Cf. swarm intelligence simple global behavior
  from complex components)
• Artificial Biochemistry
• Stochastic collectives with Law of Mass
  Interaction kinetics
• Connections to classical Markov theory, chemical
  Master Equation, and Rate Equation
• Properties of collective behavior
• Simulation
• Systematic translation to ODEs from parametric
  process libraries
• Correspondence (or not) between stochastic and
  deterministic behavior
• Interdisciplinary connections
• Process descriptions vs. chemical descriptions
• Process descriptions vs. ODE descriptions

53
Q?
54
Groupies ODEs
55
Celebrity ODE
Q What does this do?
A !a(r)A Å ?a(r)B B !b(s)B Å ?b(s)A
A B 0 hence ABN.
A -r/2(A2-B2)
A -r/2(A2(k-A)2)
AA r AB BB s BA
A -r2A(A-1) s2(N-A)(N-1-A)
-r2A2r2A s2N(N-1)-s2NA-s2(N-1)As2A2
(s2-r2)A2
A -r/2(-k22kA)
A -rk/2(-k2A)
A -rA(A-1)/2 sB(B-1)/2 B
-sB(B-1)/2 rA(A-1)/2
A -rk/2(-(AB)2A)
Stochastic Answer convergence and stability
Deterministic Answer convergence and stability
A -rk/2(A-B)
from A200.0, B0.0 intial conditions
from A200 B0 intial conditions
directive sample 0.1 1000 directive plot B()
A() new a_at_1.0chan() new b_at_1.0chan() let A()
do !a A() or ?a B() and B() do !b B() or ?b
A() run 200 of A()
interval/step 00.00030.03 r1.0 (A) dx1/dt
-x1(x1-1)/2x2(x2-1)/2 200.0 (B) dx2/dt
-x2(x2-1)/2x1(x1-1)/2 0.0
-1.0(y1-y2), 1.0 1.0(y1-y2), 0.0
Matlab continuous_sys_generator
SPiM
DESSolver
56
Groupies ODE
Same rate r
Different rates r,s
Q What does this do?
A !a(r)A Å ?b(r)B B !b(r)B Å ?a(r)A
A !a(r)A Å ?b(s)B B !b(s)B Å ?a(r)A
AB r AA BA r BB
AB r AA BA s BB
A rAB-rBA B rBA-rAB
A rAB-sBA B sBA-rAB
Stochastic Answer bistable system
A 0 B 0
A (r-s)AB B (s-r)BA
Deterministic Answer monostable systems
Deterministic Answer constant system
SPiM
ODE predicts A will always win, if rgts, and B
will always loose. The stochastic system can have
either outcome (with different probabilities).
Matlab continuous_sys_generator
ODE predicts stability A0 for any value of
A, while the stochastic system is stable only
when A is either 0 or Max.
directive sample 5.0 1000 directive plot A()
B() new a_at_1.0chan() new b_at_1.0chan() let A()
do !a A() or ?b B() and B() do !b B() or ?a
A() run 100 of (A() B())
stop
00.010.5dx1/dt 0.1x1x2, 100 dx2/dt
-0.1x1x2, 100
57
Doped Groupies ODE
Q What does this do?
directive sample 10.0 1000 directive plot Ga()
Gb() Da() Db() new a_at_1.0chan() new
b_at_1.0chan() let Ga() do !a Ga() or ?b
Gb() and Gb() do !b Gb() or ?a Ga() let Da()
!a Da() and Db() !b Db() run 1 of (Da()
Db()) run 100 of (Ga() Gb())
A !a(r)A Å ?b(r)B B !b(r)B Å ?a(r)A
Ad !a(r)Ad Bd !b(r)Bd
AB r AA BA r BB
ABd r BBd BAd r AAd
A rAB-rBA-rABdrBAd B
rBA-rAB-rBAdrABd
Ad 0 Bd 0
Ad,Bd are constant assume them both k
At B0 A-rkA, BrkA At
AB AB0 At AB AB0
A -rk(A-B) B rk(A-B)
Stochastic Answer bounded random walk
Deterministic Answer convergence and stability
-1.0(y1-y2), 1.0 1.0(y1-y2), 0.0
00.0033.0 r1.0 k1.0 dx1/dt -(x1-x2),
200.0 dx2/dt (x1-x2), 0.0
ODE predicts converging stable equilibrium at
AB instead of the total chaos observed in
the stochastic system!
For k0 (no dope), predicts deadlock AB0
but at any value of A, which is definitely not
true in the stochastic system.
Matlab continuous_sys_generator
DESSolver
SPiM
58
Hysteric Groupies ODE
Q What does this do?
directive sample 10.0 1000 directive plot Ga()
Gb() new a_at_1.0chan() new b_at_1.0chan() let Ga()
do !a Ga() or ?b ?b Gb() and Gb() do !b
Gb() or ?a ?a Ga() let Da() !a Da() and
Db() !b Db() run 100 of (Ga() Gb()) run 1
of (Da() Db())
A !a(r)A Å ?bA A ?bB B !b(r)B Å
?aB B ?aA
Ad !a(r)Ad Bd !b(r)Bd
AB r AB AB r AA BA r BA BA r
BB
ABd r ABd ABd r BBd BAd r BAd
BAd r ABd
A rAB-rBA-rABdrBAd A
rBA-rBArABd-rABd B
rBA-rAB-rBAdrABd B
rAB-rABrBAd-rBAd
A rAB-rBA-rkArkB A
rBA-rBArkA-rkA B
rBA-rAB-rkBrkA B
rAB-rABrkB-rkB
Ad 0 Bd 0
Ad,Bd are constant assume them both k
ODE predicts dampened oscillation, while the
stochasic system keeps oscillating at max level.
r1.0k1.0
Deterministic Answer dampened oscillation
Stochastic Answer robust quasi-oscillation
dx1/dtx1x4-x3x1-x1x4, 200.0 dx2/dtx3x1-x3x2
x1-x2, 0.0 dx3/dtx3x2-x1x3-x3x2,
0.0 dx4/dtx1x3-x1x4x3-x4, 0.0
SPiM
Matlab continuous_sys_generator
59
Collective Boolean Logic
60
The Strength of Populations
directive sample 0.01 1000 directive plot
B() val lam 1000.0 val mu 1.0 new
a_at_muchan let A() !a A() and B() ?a C()
and C() delay_at_lam B() run 1000 of (A()
B())
At size 2N, on a shared channel, µ is N times
stronger than ? interaction easily wins over
delay.
C
?a
!a
_at_?
_at_µ
B
A
fight!
fight!
N
N
61
Boolean Inverter Collectives
b not a
b not a
b not a
b not a
c not b
signal
no signal
in presence of a, b goes lowin absence of a, b
goes high
the high b state reinforces itself (as a
population)
input stimulus
perfect rectifier
hysteresis
zero-point noise resistant
directive sample 110.0 1000 directive plot !a
!b new a_at_1.0chan new b_at_1.0chan let
Inv2_hi(achan, bchan) do !b Inv2_hi(a,b)
or ?a Inv2_mi(a,b) and Inv2_mi(achan, bchan)
do ?b Inv2_hi(a,b) or delay_at_1.0
Inv2_hi(a,b) or ?a Inv2_lo(a,b) and
Inv2_lo(achan, bchan) do ?b Inv2_mi(a,b)
or delay_at_1.0 Inv2_mi(a,b) run 100 of
Inv2_hi(a,b) let clock(tfloat, tickchan)
( sends a tick every t time ) (val ti
t/100.0 val d 1.0/ti ( by 100-step erlang
timers ) let step(nint) if nlt0 then
!tick clock(t,tick) else delay_at_d step(n-1)
run step(100)) let S1(achan, tockchan) do
!a S1(a,tock) or ?tock () let SN(nint,
tfloat, achan, tickchan, tockchan) if
n0 then clock(t, tock) else ?tick (S1(a,tock)
SN(n-1,t,a,tick,tock)) let raisingfalling(achan,
nint, tfloat) (new tickchan new tockchan
run (clock(t,tick) SN(n,t,a,tick,tock))) ru
n raisingfalling(a,100,0.5)
directive sample 110.0 1000 directive plot !a
!b !c !d new a_at_1.0chan new b_at_1.0chan new
c_at_1.0chan let Inv2_hi(achan, bchan) do
!b Inv2_hi(a,b) or ?a Inv2_mi(a,b) and
Inv2_mi(achan, bchan) do ?b Inv2_hi(a,b)
or delay_at_1.0 Inv2_hi(a,b) or ?a
Inv2_lo(a,b) and Inv2_lo(achan, bchan) do
?b Inv2_mi(a,b) or delay_at_1.0 Inv2_mi(a,b) run
100 of (Inv2_hi(a,b) Inv2_lo(b,c)) let
clock(tfloat, tickchan) ( sends a
tick every t time ) (val ti t/100.0 val d
1.0/ti ( by 100-step erlang timers ) let
step(nint) if nlt0 then !tick clock(t,tick)
else delay_at_d step(n-1) run step(100)) let
S1(achan, tockchan) do !a S1(a,tock) or
?tock () let SN(nint, tfloat, achan,
tickchan, tockchan) if n0 then clock(t,
tock) else ?tick (S1(a,tock)
SN(n-1,t,a,tick,tock)) let raisingfalling(achan,
nint, tfloat) (new tickchan new tockchan
run (clock(t,tick) SN(n,t,a,tick,tock))) ru
n raisingfalling(a,100,0.5)
directive sample 110.0 1000 directive plot !a
!b new a_at_1.0chan new b_at_1.0chan let
Inv_hi(achan, bchan) do !b Inv_hi(a,b)
or ?a Inv_lo(a,b) and Inv_lo(achan, bchan)
do ?b Inv_hi(a,b) or delay_at_1.0
Inv_hi(a,b) run 100 of Inv_hi(a,b) let
clock(tfloat, tickchan) ( sends a
tick every t time ) (val ti t/100.0 val d
1.0/ti ( by 100-step erlang timers ) let
step(nint) if nlt0 then !tick clock(t,tick)
else delay_at_d step(n-1) run step(100)) let
S1(achan, tockchan) do !a S1(a,tock) or
?tock () let SN(nint, tfloat, achan,
tickchan, tockchan) if n0 then clock(t,
tock) else ?tick (S1(a,tock)
SN(n-1,t,a,tick,tock)) let raisingfalling(achan,
nint, tfloat) (new tickchan new tockchan
run (clock(t,tick) SN(n,t,a,tick,tock))) ru
n raisingfalling(a,100,0.5)
directive sample 110.0 1000 directive plot !a
!b new a_at_1.0chan new b_at_1.0chan let
Inv_hi(achan, bchan) do !b Inv_hi(a,b)
or ?a Inv_lo(a,b) and Inv_lo(achan, bchan)
delay_at_1.0 Inv_hi(a,b) run 100 of
Inv_hi(a,b) let clock(tfloat, tickchan)
( sends a tick every t time ) (val ti
t/100.0 val d 1.0/ti ( by 100-step erlang
timers ) let step(nint) if nlt0 then
!tick clock(t,tick) else delay_at_d step(n-1)
run step(100)) let S1(achan, tockchan) do
!a S1(a,tock) or ?tock () let SN(nint,
tfloat, achan, tickchan, tockchan) if
n0 then clock(t, tock) else ?tick (S1(a,tock)
SN(n-1,t,a,tick,tock)) let raisingfalling(achan,
nint, tfloat) (new tickchan new tockchan
run (clock(t,tick) SN(n,t,a,tick,tock))) ru
n raisingfalling(a,100,0.5)
62
Boolean Gate Collectives
c a or b
c a imply b
c a unless b
c a and b
c a xor b
Inputs 10 !a for 4t 2t 10 !b for 4t
directive sample 10.0 1000 directive plot !a !b
!c new a_at_1.0chan new b_at_1.0chan new
c_at_1.0chan let Xor_hi_a(achan, bchan, cchan)
do !c Xor_hi_a(a,b,c) or ?b
Xor_lo_ab(a,b,c) or delay_at_1.0 Xor_lo_a(a,b,c)
and Xor_hi_b(achan, bchan, cchan) do !c
Xor_hi_b(a,b,c) or ?a Xor_lo_ab(a,b,c) or
delay_at_1.0 Xor_lo_b(a,b,c) and Xor_lo_a(achan,
bchan, cchan) do ?a Xor_hi_a(a,b,c) or
?b Xor_lo_ab(a,b,c) and Xor_lo_b(achan, bchan,
cchan) do ?b Xor_hi_b(a,b,c) or ?a
Xor_lo_ab(a,b,c) and Xor_lo_ab(achan, bchan,
cchan) do delay_at_1.0 Xor_hi_a(a,b,c) or
delay_at_1.0 Xor_hi_b(a,b,c) run 50 of
(Xor_lo_a(a,b,c) Xor_lo_b(a,b,c)) let
clock(tfloat, tickchan) ( sends a tick
every t time ) (val ti t/200.0 val d
1.0/ti let step(nint) if nlt0 then
!tick clock(t, tick) else delay_at_d step(n-1)
run step(200)) let S_a(tickchan) do !a
S_a(tick) or ?tick () let S_b(tickchan)
?tick S_b1(tick) and S_b1(tickchan) do !b
S_b1(tick) or ?tick S_b2(tick) and
S_b2(tickchan) do !b S_b2(tick) or ?tick
() run 10 of (new tickchan run (clock(4.0,tick)
S_a(tick))) run 10 of (new tickchan run
(clock(2.0,tick) S_b(tick)))
directive sample 10.0 1000 directive plot !a !b
!c new a_at_1.0chan new b_at_1.0chan new
c_at_1.0chan val del 1.0 let And_hi(achan,
bchan, cchan) do !c And_hi(a,b,c) or
delay_at_del And_lo_a(a,b,c) and And_lo_a(achan,
bchan, cchan) do ?a And_hi(a,b,c) or
delay_at_del And_lo_b(a,b,c) and And_lo_b(achan,
bchan, cchan) ?b And_lo_a(a,b,c) run 100
of And_lo_b(a,b,c) let clock(tfloat, tickchan)
( sends a tick every t time ) (val ti
t/200.0 val d 1.0/ti let step(nint)
if nlt0 then !tick clock(t, tick) else
delay_at_d step(n-1) run step(200)) let
S_a(tickchan) do !a S_a(tick) or ?tick
() let S_b(tickchan) ?tick S_b1(tick) and
S_b1(tickchan) do !b S_b1(tick) or ?tick
S_b2(tick) and S_b2(tickchan) do !b
S_b2(tick) or ?tick S_b3(tick) and
S_b3(tickchan) do !b S_b3(tick) or ?tick
() run 10 of (new tickchan run (clock(4.0,tick)
S_a(tick))) run 10 of (new tickchan run
(clock(2.0,tick) S_b(tick)))
directive sample 10.0 1000 directive plot !a !b
!c new a_at_1.0chan new b_at_1.0chan new
c_at_1.0chan val del 1.0 let OOlO_hi(achan,
bchan, cchan) do !c OOlO_hi(a,b,c) or
delay_at_del OOlO_lo_a(a,b,c) or ?b
OOlO_lo_b(a,b,c) and OOlO_lo_a(achan, bchan,
cchan) ?a OOlO_hi(a,b,c) and
OOlO_lo_b(achan, bchan, cchan) delay_at_del
OOlO_hi(a,b,c) run 50 of (OOlO_lo_a(a,b,c)
OOlO_lo_b(a,b,c)) let clock(tfloat, tickchan)
( sends a tick every t time ) (val ti
t/200.0 val d 1.0/ti let step(nint)
if nlt0 then !tick clock(t, tick) else
delay_at_d step(n-1) run step(200)) let
S_a(tickchan) do !a S_a(tick) or ?tick
() let S_b(tickchan) ?tick S_b1(tick) and
S_b1(tickchan) do !b S_b1(tick) or ?tick
S_b2(tick) and S_b2(tickchan) do !b
S_b2(tick) or ?tick () run 10 of (new tickchan
run (clock(4.0,tick) S_a(tick))) run 10 of (new
tickchan run (clock(2.0,tick) S_b(tick)))
directive sample 10.0 1000 directive plot !a !b
!c new a_at_1.0chan new b_at_1.0chan new
c_at_1.0chan val del 1.0 let Imply_hi_a(achan,
bchan, cchan) do !c Imply_hi_a(a,b,c) or
?a Imply_lo(a,b,c) and Imply_hi_b(achan,
bchan, cchan) do !c Imply_hi_b(a,b,c) or
delay_at_del Imply_lo(a,b,c) and Imply_lo(achan,
bchan, cchan) do ?b Imply_hi_b(a,b,c) or
delay_at_del Imply_hi_a(a,b,c) run 100 of
Imply_lo(a,b,c) let clock(tfloat, tickchan)
( sends a tick every t time ) (val ti
t/200.0 val d 1.0/ti let step(nint)
if nlt0 then !tick clock(t, tick) else delay_at_d
step(n-1) run step(200)) let S_a(tickchan)
do !a S_a(tick) or ?tick () let S_b(tickchan)
?tick S_b1(tick) and S_b1(tickchan) do !b
S_b1(tick) or ?tick S_b2(tick) and
S_b2(tickchan) do !b S_b2(tick) or ?tick
() run 10 of (new tickchan run (clock(4.0,tick)
S_a(tick))) run 10 of (new tickchan run
(clock(2.0,tick) S_b(tick)))
directive sample 10.0 1000 directive plot !a !b
!c new a_at_1.0chan new b_at_1.0chan new
c_at_1.0chan val del 1.0 let Or_hi(achan,
bchan, cchan) do !c Or_hi(a,b,c) or
delay_at_del Or_lo(a,b,c) and Or_lo(achan, bchan,
cchan) do ?a Or_hi(a,b,c) or ?b
Or_hi(a,b,c) run 100 of Or_lo(a,b,c) let
clock(tfloat, tickchan) ( sends a tick
every t time ) (val ti t/200.0 val d
1.0/ti let step(nint) if nlt0 then
!tick clock(t, tick) else delay_at_d step(n-1)
run step(200)) let S_a(tickchan) do !a
S_a(tick) or ?tick () let S_b(tickchan)
?tick S_b1(tick) and S_b1(tickchan) do !b
S_b1(tick) or ?tick S_b2(tick) and
S_b2(tickchan) do !b S_b2(tick) or ?tick
() run 10 of (new tickchan run (clock(4.0,tick)
S_a(tick))) run 10 of (new tickchan run
(clock(2.0,tick) S_b(tick)))
63
Collective Analog Devices
64
Xor as an Op Amp
c A(a b) d A(b - a)
Follower (a standard OpAmp trick)
(!c!a)
?a
a0 b0 Þ db-a0 aca-b0 a0 b1 Þ db-a1
aca-b0 a1 b0 Þ db-a0 aca-b1 a1 b1 Þ
db-a0 aca-b0 hence d1 at next step
!d
?b
Noninverting Configuration
hence db
db analog response!!
directive sample 20.0 1000 directive plot !a !b
!c !d new a_at_1.0chan new b_at_1.0chan new
c_at_1.0chan new d_at_1.0chan let Xor_hi_a(achan,
bchan, cchan, dchan) do !c
Xor_hi_a(a,b,c,d) or ?b Xor_lo_ab(a,b,c,d) or
delay_at_1.0 Xor_lo_a(a,b,c,d) and
Xor_hi_b(achan, bchan, cchan, dchan) do
!d Xor_hi_b(a,b,c,d) or ?a Xor_lo_ab(a,b,c,d)
or delay_at_1.0 Xor_lo_b(a,b,c,d) and
Xor_lo_a(achan, bchan, cchan, dchan) do
?a Xor_hi_a(a,b,c,d) or ?b Xor_lo_ab(a,b,c,d) an
d Xor_lo_b(achan, bchan, cchan, dchan)
do ?b Xor_hi_b(a,b,c,d) or ?a
Xor_lo_ab(a,b,c,d) and Xor_lo_ab(achan, bchan,
cchan, dchan) do delay_at_1.0
Xor_hi_a(a,b,c,d) or delay_at_1.0
Xor_hi_b(a,b,c,d) run 50 of (Xor_lo_a(a,b,c,d)
Xor_lo_b(a,b,c,d)) let clock(tfloat, tickchan)
( sends a tick every t time ) (val ti
t/200.0 val d 1.0/ti let step(nint)
if nlt0 then !tick clock(t, tick) else
delay_at_d step(n-1) run step(200)) let
S_a(tickchan) do !a S_a(tick) or ?tick
() let S_b(tickchan) ?tick S_b1(tick) and
S_b1(tickchan) do !b S_b1(tick) or ?tick
S_b2(tick) and S_b2(tickchan) do !b
S_b2(tick) or ?tick () run 100 of (new
tickchan run (clock(8.0,tick) S_a(tick))) run
100 of (new tickchan run (clock(4.0,tick)
S_b(tick)))
directive sample 40.0 1000 directive plot !a !b
!d new a_at_1.0chan new b_at_1.0chan new
d_at_1.0chan let Xor_hi_a(achan, bchan, cchan,
dchan) do !c Xor_hi_a(a,b,c,d) or ?b
Xor_lo_ab(a,b,c,d) or delay_at_1.0
Xor_lo_a(a,b,c,d) and Xor_hi_b(achan, bchan,
cchan, dchan) do !d Xor_hi_b(a,b,c,d) or
?a Xor_lo_ab(a,b,c,d) or delay_at_1.0
Xor_lo_b(a,b,c,d) and Xor_lo_a(achan, bchan,
cchan, dchan) do ?a Xor_hi_a(a,b,c,d) or
?b Xor_lo_ab(a,b,c,d) and Xor_lo_b(achan,
bchan, cchan, dchan) do ?b
Xor_hi_b(a,b,c,d) or ?a Xor_lo_ab(a,b,c,d) and
Xor_lo_ab(achan, bchan, cchan, dchan) do
delay_at_1.0 Xor_hi_a(a,b,c,d) or delay_at_1.0
Xor_hi_b(a,b,c,d) run 50 of (Xor_lo_a(a,b,a,d)
Xor_lo_b(a,b,a,d)) let clock(tfloat, tickchan)
( sends a tick every t time ) (val ti
t/200.0 val d 1.0/ti let step(nint)
if nlt0 then !tick clock(t, tick) else
delay_at_d step(n-1) run step(200)) let
S_b(tickchan) ?tick S_b1(tick) and
S_b1(tickchan) do !b S_b1(tick) or ?tick
S_b2(tick) and S_b2(tickchan) do !b
S_b2(tick) or ?tick S_b3(tick) and
S_b3(tickchan) ?tick S_b4(tick) and
S_b4(tickchan) !b S_b4(tick) run 10 of (new
tickchan run (clock(8.0,tick) S_b(tick)))
a100 may or may not happen
65
Changing the OpAmp Gain
An OpAmp provides infinite differential
amplification, but a stable finite amplification
can be obtained by a feedback loop with a load
splitter (the follower is a special case of that,
which gives gain 1). The equivalent here is
simply changing the rate on the feedback link.
rate control
?a
!d
?b
Empirical lawd b/rate(a) but why?
directive sample 40.0 1000 directive plot !a !b
!d new a_at_1.0chan new b_at_1.0chan new
d_at_1.0chan let Xor_hi_a(achan, bchan, cchan,
dchan) do !c Xor_hi_a(a,b,c,d) or ?b
Xor_lo_ab(a,b,c,d) or delay_at_1.0
Xor_lo_a(a,b,c,d) and Xor_hi_b(achan, bchan,
cchan, dchan) do !d Xor_hi_b(a,b,c,d) or
?a Xor_lo_ab(a,b,c,d) or delay_at_1.0
Xor_lo_b(a,b,c,d) and Xor_lo_a(achan, bchan,
cchan, dchan) do ?a Xor_hi_a(a,b,c,d) or
?b Xor_lo_ab(a,b,c,d) and Xor_lo_b(achan,
bchan, cchan, dchan) do ?b
Xor_hi_b(a,b,c,d) or ?a Xor_lo_ab(a,b,c,d) and
Xor_lo_ab(achan, bchan, cchan, dchan) do
delay_at_1.0 Xor_hi_a(a,b,c,d) or delay_at_1.0
Xor_hi_b(a,b,c,d) run 100 of (Xor_lo_a(a,b,a,d)
Xor_lo_b(a,b,a,d)) run 100 of replicate !b
b100 a_at_0.5 d gain 2.00OpAmp400 (non saturated)
b100 a_at_1.0 d gain 1.0 OpAmp200
b100 a_at_0.6 d gain 1.66OpAmp200
b100 a_at_0.75 d gain 1.33OpAmp200
b100 a_at_2.0 d gain 0.5OpAmp200
b100 a_at_0.5 d gain 2.00OpAmp200 (saturated)
b100 a_at_4.0 d gain 0.25OpAmp200
b100 a_at_0.33 d gain 3.00OpAmp400
66
Automata Polymers
67
BidirectionalPolymerization
A
A
A
A
new c_at_µ new stop_at_1.0 Afree (new rht_at_?
!c(rht) Abrht(rht)) ?c(lft)
Ablft(lft) Ablft(lft) (new rht_at_? !c(rht)
Abound(lft,rht)) Abrht(rht) ?c(lft)
Abound(lft,rht) Abound(lft,rht) ?stop
?c(r)
!c(nl)
Monomer Automaton
Communicating Automata Bound output !c(nr) and
input ?c(l) on automata transitions to model
complexation
Afree
Ablft
Abrht
Abound
directive sample 10000.0 directive plot Afree()
Ablft() Abrht() Abound() val lam 1.0 val
mu 1.0 new c_at_muchan(chan) new
stop_at_1.0chan let Afree() (new rht_at_lamchan
run do !c(rht) Abrht(rht) or ?c(lft)
Ablft(lft)) and Ablft(lftchan) (new
rht_at_lamchan run !c(rht) Abound(lft,rht)) and
Abrht(rhtchan) ?c(lft)
Abound(lft,rht) and Abound(lftchan, rhtchan)
?stop run (2 of Afree())
68
Bidirectional PolymerizationCircular Polymer
Lengths
Scanning and counting the size of the circular
polymers (by a cheap trick). Polymer formation is
complete within 10t then a different polymer is
scanned every 100t.
directive sample 1000.0 directive plot Abound()
?count type Link chan(chan) type Barb
chan val lam 1000.0 ( set high for better
counting ) val mu 1.0 new c_at_muchan(Link) new
enter_at_lamchan(Barb) new count_at_lamBarb let
Afree() (new rht_at_lamLink run do !c(rht)
Abrht(rht) or ?c(lft) Ablft(lft)) and
Ablft(lftLink) (new rht_at_lamLink run
!c(rht) Abound(lft,rht)) and Abrht(rhtLink)
?c(lft) Abound(lft,rht) and
Abound(lftLink, rhtLink) do ?enter(barb)
(?barb !rht(barb)) or ?lft(barb) (?barb
!rht(barb)) ( each Abound waits for a barb,
exhibits it, and passes it to the right so we
can plot number of Abound in a ring ) let
clock(tfloat, tickchan) ( sends a tick
every t time ) (val ti t/1000.0 val d
1.0/ti let step(nint) if nlt0 then
!tick clock(t,tick) else delay_at_d step(n-1)
run step(1000)) new tickchan let Scan()
?tick !enter(count) Scan() run 100 of Afree()
run (clock(100.0, tick) Scan())
100xAfree, initially. The height of each rising
step is the size of a separate circular polymer.
(Unbiased sample of nine consecutive runs.)
69
Actin-likePoly/Depolymerization
new c_at_µ Afree (new lft_at_? !c(lft)
Ablft(lft)) ?c(rht) Abrht(rht) Ablft(lft)
!lft Afree ?c(rht) Abound(lft,rht) Abr
ht(rht) ?rht Afree Abound(lft,rht) !lft
Abrht(rht)
!c(nl)
?c(r)
?r
!l
Monomer Automaton
Free
?c(r)
!c(nl)
?c(r)
!c(nl)
rnl
!l
?r
!l
?r
Bound left
Boundright
!l
!l
Bound right
Bound left
?c(r)
?c(r)
Boundboth
Boundboth
70
(Note on computing the multisets)

• A multiset M Î Multiset(S), where S is a set with
  equality, is a total function SNat, which may
  also be written as a finite enumeration with
  repetitions m.
• Multiset binary union is the function Èm(M,M)
  \s. M(s)M(s).
• The shorthand
• (X r P) s.t. (X º trP Å ) Î Em
• is defined as the following finite union of
  singleton multisets
• Èm(Xp1P1ÅÅpnPn)ÎE
• of(Èmi s.t. pitr
• of (X r Pi)m)
• i.e. for each (Xp1P1ÅÅpnPn)ÎE and for each i
  such that pitr, return a copy of (X r Pi).
• The shorthand
• (X Y r P Q) s.t. X?Y, á(X º ?n(r)P Å
  ),(Y º !n(r)Q Å )ñ Î E2m

71
Possible Interactions
The speed of interaction is proportional to the
number of possible interactions. But a process
cannot interact with itself.
2
?2a
!1a
In(a,A) 2 Out(a,A) 3 Mix(a,A) 23 6
In(a,B) 2 Out(a,B) 1 Mix(a,B) 21 2
Assume each process P is in restricted-sum-normal-
form. For each channel x In(x,P) Num of
active ?x in P Out(x,P) Num of active !x in
P Mix(x,P) In(x,P)Out(x,P) In(x) Sum P of
In(x,P) Out(x) Sum P of Out(x,P) Mix(x) Sum P
of Mix(x,P) The global Activity on channel
x The global speed of interaction on a
channel x
B
6
2
!3a
?2a
6
Forbidden (Mix) interactions
I