Computational Methods in Systems Biology and Synthetic Biology Franois Fages, Constraint Programming

1
Computational Methods inSystems Biology and
Synthetic BiologyFrançois Fages, Constraint
Programming Group, INRIA Rocquencourt
mailtoFrancois.Fages_at_inria.frhttp//contraintes.
inria.fr/
2
Overview of the Lectures

• Formal molecules and reaction models in BIOCHAM
• Kinetics
• Qualitative properties formalized in temporal
  logic CTL
• Quantitative properties formalized in LTL(R) and
  pLTL(R)

3
Biochemical Kinetics

• Study the concentration of chemical substances in
  a biological system as a function of time.
• Continuous semantics of Biocham
• Molecules A1 ,, Am
  
• ANumber of molecules A
• AConcentration of A in the solution A A
  / Volume ML-1
• Solutions with stoichiometric coefficients c1
  A1 cn An

4
Law of Mass Action

• Assumption each molecule moves independently of
  other molecules in a random walk (diffusion,
  dilute solutions, low concentration)
• The number of interactions of A with B is
  proportional to the number of A and B molecules,
  the proportionality factor k is the rate constant
  of the reaction
• A B ?k C
• the rate of the reaction is kAB.
• dC/dt k A B
• dA/dt .

5
Law of Mass Action

• Assumption each molecule moves independently of
  other molecules in a random walk (diffusion,
  dilute solutions, low concentration)
• The number of interactions of A with B is
  proportional to the number of A and B molecules,
  the proportionality factor k is the rate constant
  of the reaction
• A B ?k C
• the rate of the reaction is kAB
• dC/dt k A B
• dA/dt -k A B
• dB/dt -k A B

6
Interpretation of Rate Constants k

• Complexation probability of reaction upon
  collision (specificity, affinity)
• position of matching surfaces
• Decomplexation total energy of all bonds
• (giving dissociation rates)
• Different diffusion speeds (small
  moleculesgtsubstratesgtenzymes)
• Average travel in a random walk 1 µm in 1s,
  2µm in 4s, 10µm in 100s
• For one enzyme
• 500000 random collisions per second with
  substrate concentration of 10-5
• 50000 random collisions per second with
  substrate concentration of 10-6

7
BIOCHAM Concentration Semantics

• To a set of BIOCHAM rules with kinetic
  expressions ei
• ei for SigtSI i1,,n
• one associates the system of ODEs over variables
  A1 ,, Ak
• dAk/dt Sni1 ri(Ak)ei - Snj1 lj(Ak)ej
• where ri(A) is the stoichiometric coefficient
  of A in Si
• li(A) is the stoichiometric
  coefficient of A in Si .

8
Signal Reception on the Membrane

present(L,0.5). present(RTK,0.01). absent(L-RTK).
absent(S). parameter(k1,1). parameter(k2,0.1). pa
rameter(k3,1). parameter(k4,0.3). (k1LRTK,
k2L-RTK) for LRTK ltgt L-RTK. (k3L-RTK2,
k4S) for 2(L-RTK) ltgt S. dS/dT
9
Signal Reception on the Membrane

present(L,0.5). present(RTK,0.01). absent(L-RTK).
absent(S). parameter(k1,1). parameter(k2,0.1). pa
rameter(k3,1). parameter(k4,0.3). (k1LRTK,
k2L-RTK) for LRTK ltgt L-RTK. (k3L-RTK2,
k4S) for 2(L-RTK) ltgt S. dS/dT
k3L-RTK2 k4S d(L-RTK)/dT
10
Signal Reception on the Membrane

present(L,0.5). present(RTK,0.01). absent(L-RTK).
absent(S). parameter(k1,1). parameter(k2,0.1). pa
rameter(k3,1). parameter(k4,0.3). (k1LRTK,
k2L-RTK) for LRTK ltgt L-RTK. (k3L-RTK2,
k4S) for 2(L-RTK) ltgt S. dS/dT
k3L-RTK2 k4S d(L-RTK)/dT k1LRTK
2k4S k2L-RTK 2k3L-RTK2
11
Compositionality of Reaction Rules

• The union of two sets of reaction rules is a set
  of reaction rules
• Rule-based models can thus be composed to form
  complex reaction models by set union
  (alternative addition of the kinetics)
• Modular models (e.g. for synthetic biology) would
  need

12
Compositionality of Reaction Rules

• The union of two sets of reaction rules is a set
  of reaction rules
• Rule-based models can thus be composed to form
  complex reaction models by set union
  (alternative addition of the kinetics)
• Modular models (e.g. for synthetic biology) would
  need
• Sufficiently decomposed reaction rules
• ESltgtC gtEP. (not SltEgtP if
  competition on C)
• Sufficiently general kinetics expression
• parameters as possibly functions of temperature,
  pH, pressure, light,
• different pH-logH in intracellular and
  extracellular solvents (water)
• Ex. pH(cytosol)7.2, pH(lysosomes)4.5,
  pH(cytoplasm) in 6.6,7.2
• Interface variables
• controlled by other modules (endogeneous
  variables)
• controlled by external laws (exogeneous
  variables).

13
Michaelis-Menten Enzymatic Reaction

• An enzyme E binds to a substrate S to catalyze
  the formation of product P
• ES ?k1 C ?k2 EP
• ES ?km1 C
• Compiles into a system of non-linear Ordinary
  Differential Equations
• dE/dt -k1ES(k2km1)C
• dS/dt

14
Michaelis-Menten Enzymatic Reaction

• An enzyme E binds to a substrate S to catalyze
  the formation of product P
• ES ?k1 C ?k2 EP
• ES ?km1 C
• Compiles into a system of non-linear Ordinary
  Differential Equations
• dE/dt -k1ES(k2km1)C
• dS/dt -k1ESkm1C
• dC/dt

15
Michaelis-Menten Enzymatic Reaction

• An enzyme E binds to a substrate S to catalyze
  the formation of product P
• ES ?k1 C ?k2 EP
• ES ?km1 C
• Compiles into a system of non-linear Ordinary
  Differential Equations
• dE/dt -k1ES(k2km1)C
• dS/dt -k1ESkm1C
• dC/dt k1ES-(k2km1)C
• dP/dt k2C

16
Michaelis-Menten Enzymatic Reaction

• An enzyme E binds to a substrate S to catalyze
  the formation of product P
• ES ?k1 C ?k2 EP
• ES ?km1 C
• Compiles into a system of non-linear Ordinary
  Differential Equations
• dE/dt -k1ES(k2km1)C
• dS/dt -k1ESkm1C
• dC/dt k1ES-(k2km1)C
• dP/dt k2C
• Assuming C0P00, we get EE0-C, S0SCP,
• dS/dt -k1(E0-C)Skm1C
• dC/dt k1(E0-C)S-(k2km1)C

17
Numerical Integration Methods

• System dX/dt f(X).
• Initial conditions X0
• Idea discretize time t0, t1t0?t, t2t1?t,
• and compute a trace
• (t0,X0,dX0/dt), (t1,X1,dX1/dt), ,
  (tn,Xn,dXn/dt)

18
Numerical Integration Methods

• System dX/dt f(X).
• Initial conditions X0
• Idea discretize time t0, t1t0?t, t2t1?t,
• and compute a trace
• (t0,X0,dX0/dt), (t1,X1,dX1/dt), ,
  (tn,Xn,dXn/dt)
• Eulers method ti1ti ?t
• Xi1Xif(Xi)?t
• error estimation E(Xi1)f(Xi)-f(Xi1)?t

19
Numerical Integration Methods

• System dX/dt f(X).
• Initial conditions X0
• Idea discretize time t0, t1t0?t, t2t1?t,
• and compute a trace
• (t0,X0,dX0/dt), (t1,X1,dX1/dt), ,
  (tn,Xn,dXn/dt)
• Eulers method ti1ti ?t
• Xi1Xif(Xi)?t
• error estimation E(Xi1)f(Xi)-f(Xi1)?t
• Runge-Kuttas method intermediate computations
  at ?t/2

20
Numerical Integration Methods

• System dX/dt f(X).
• Initial conditions X0
• Idea discretize time t0, t1t0?t, t2t1?t,
• and compute a trace
• (t0,X0,dX0/dt), (t1,X1,dX1/dt), ,
  (tn,Xn,dXn/dt)
• Eulers method ti1ti ?t
• Xi1Xif(Xi)?t
• error estimation E(Xi1)f(Xi)-f(Xi1)?t
• Runge-Kuttas method intermediate computations
  at ?t/2
• Adaptive step method ?ti1 ?ti/2 while EgtEmax,
  otherwise ?ti1 2?ti

21
Numerical Integration Methods

• System dX/dt f(X).
• Initial conditions X0
• Idea discretize time t0, t1t0?t, t2t1?t,
• and compute a trace
• (t0,X0,dX0/dt), (t1,X1,dX1/dt), ,
  (tn,Xn,dXn/dt)
• Eulers method ti1ti ?t
• Xi1Xif(Xi)?t
• error estimation E(Xi1)f(Xi)-f(Xi1)?t
• Runge-Kuttas method intermediate computations
  at ?t/2
• Adaptive step method ?ti1 ?ti/2 while EgtEmax,
  otherwise ?ti1 2?ti
• Rosenbrocks stiff method solve
  Xi1Xif(Xi1)?t by formal differentiation

22
Multi-Scale Phenomena
Hydrolysis of benzoyl-L-arginine ethyl ester by
trypsin present(En,1e-8). present(S,1e-5).
absent(C). absent(P). parameter(k1,4e6).
parameter(km1,25). parameter(k2,15). (k1EnS,
km1C) for EnS ltgt C. k2C for C
gt EnP. Complex formation 5e-9 in 0.1s
Product formation 1e-5 in 1000s
23
Quasi-Steady State Approximation

• After short initial period (0.1s) the complex
  concentration reaches its limit.
• Assume dC/dt0

24
Quasi-Steady State Approximation

• After short initial period (0.1s) the complex
  concentration reaches its limit.
• Assume dC/dt0
• From dC/dt k1S(E0-C)-(k2km1)C
• we get C k1E0S/(k2km1k1S)

25
Quasi-Steady State Approximation

• After short initial period (0.1s) the complex
  concentration reaches its limit.
• Assume dC/dt0
• From dC/dt k1S(E0-C)-(k2km1)C
• we get C k1E0S/(k2km1k1S)
• E0S/(((k2km1)/k1)S)
• E0S/(KmS)
• where Km(k2km1)/k1

26
Quasi-Steady State Approximation

• After short initial period (0,1s) the complex
  concentration reaches its limit.
• Assume dC/dt0
• From dC/dt k1S(E0-C)-(k2km1)C
• we get C k1E0S/(k2km1k1S)
• E0S/(((k2km1)/k1)S)
• E0S/(KmS)
• where Km(k2km1)/k1
• dS/dt -dP/dt
• -k2C
• -VmS / (KmS)
• where Vm k2E0.

27
Quasi-Steady State Approximation

• Assuming dC/dt0, we have dE/dt0 and C E0 S /
  (KmS).
• Michaelis-Menten rate dP/dt -dS/dt VmS /
  (KmS) is reaction velocity
• Vmk2E0
• Km(km1k2)/k1

28
Quasi-Steady State Approximation

• Assuming dC/dt0, we have dE/dt0 and C E0 S /
  (KmS).
• Michaelis-Menten rate dP/dt -dS/dt VmS /
  (KmS) is reaction velocity
• Vmk2E0 is maximum velocity at
  saturating substrate concentration
• Km(km1k2)/k1

29
Quasi-Steady State Approximation

• Assuming dC/dt0, we have dE/dt0 and C E0 S /
  (KmS).
• Michaelis-Menten rate dP/dt -dS/dt VmS /
  (KmS) is reaction velocity
• Vmk2E0 is maximum velocity at
  saturing substrate concentration
• Km(km1k2)/k1 is substrate concentration with
  half maximum velocity
• Experimental measurement
• The initial velocity is
• linear in E0
• hyperbolic in S0

30
Quasi-Steady State Approximation

• Assuming dC/dt0, hence dE/dt0 and C E0 S /
  (KmS).
• Michaelis-Menten rate dP/dt -dS/dt VmS /
  (KmS)
• Vmk2E0
• Km(km1k2)/k1
• macro(Vm, k2En).
• macro(Km, (km1k2)/k1).
• MM(Vm,Km) for S Engt P.
• Is equivalent to
• macro(Kf, VmS/(KmS)).
• Kf for S Engt P.

31
Competitive Inhibition

• present(En,1e-8). present(S,1e-5).
  present(I,1e-5). parameter(k3,5e5).
• parameter(k1,4e6). parameter(km1,25).
  parameter(k2,15).
• (k1EnS,km1C) for EnS ltgt C.
  k2C for C gt EnP.
• k3CI for CI gt CI.
• Complex formation 4e-9 in 0.04s
  Product formation 3e-8 in 3s

32
Isosteric Inhibition

• present(En,1e-8). present(S,1e-5).
  present(I,1e-5). parameter(k3,5e5).
• parameter(k1,4e6). parameter(km1,25).
  parameter(k2,15).
• (k1EnS,km1C) for EnS ltgt C.
  k2C for C gt EnP.
• k3EnI for EnI gt EI.
• Complex formation 2.5e-9 in 0.4s
  Product formation 2.5e-9 in 1000s

33
Allosteric Inhibition or Activation

• (iEnI,imEI) for EnI ltgt EI.
  parameter(i,1e7). parameter(im,10).
• (i1EIS,im1CI) for EIS ltgt CI.
  parameter(i1,5e6). parameter(im1,5).
• i2CI for CI gt EIP.
  parameter(i2,2).
• Complex formation 2e-9 in 0.4s Product
  formation 1e-5 in 1000s

34
Cooperative Enzymes and Hill Equation

• Dimer enzyme with two promoters
• ES ?2k1 C1 ?k2 EP C1S ?k1 C2 ?2k2
  C1P
• ES ?k-1 C1 C1S ?2k-1 C2
• Let Km(k-1k2)/k1 and Km(k-1k2)/k1
• Non-cooperative if Km Km
• Michaelis-Menten rate
  VmS / (KmS) where Vm2k2E0.
• hyperbolic velocity vs
  substrate concentration
• Cooperative if k1gtk1
• Hill equation rate
  VmS2 / (KmKmS2) where Vm2k2E0.
• sigmoid velocity vs
  substrate concentration

35
Gene Regulatory Networks Implementation

• Gene a, product Pa, promotion factor PFa (same
  for gene b)
• a PFa gt a-PFa. b PFb gt b-PFb.
• _ a-PFagt Pa. _ b-PFbgt Pb.

36
Gene Regulatory Networks Implementation

• Gene a, product Pa, promotion factor PFa (same
  for gene b)
• a PFa gt a-PFa. b PFb gt b-PFb.
• _ a-PFagt Pa. _ b-PFbgt Pb.
• a activates b a ? b (the product of a is the
  promotion factor of b)
• Pa PFb

37
Gene Regulatory Networks Implementation

• Gene a, product Pa, promotion factor PFa (same
  for gene b)
• a PFa gt a-PFa. b PFb gt b-PFb.
• _ a-PFagt Pa. _ b-PFbgt Pb.
• a activates b a ? b (the product of a is the
  promotion factor of b)
• Pa PFb
• a inhibits b a -- b (first implementation by
  transcriptional inhibition)
• b Pa gt b-Pa.

38
Gene Regulatory Networks Implementation

• Gene a, product Pa, promotion factor PFa (same
  for gene b)
• a PFa ltgt a-PFa. b PFb ltgt b-PFb.
• _ a-PFagt Pa. _ b-PFbgt Pb.
• a activates b a ? b (the product of a is the
  promotion factor of b)
• Pa PFb
• a inhibits b a -- b (first implementation by
  transcriptional inhibition)
• b Pa ltgt b-Pa.
• a inhibits b a -- b (second implementation by
  promotion factor inhibition)
• PFb Pa ltgt PFb-Pa.

39
Gene Circuit

• a activates a, a activates b, b inhibits a
• First implementation
• a Pa ltgt a-Pa.
• _ a-Pagt Pa.
• b Pa ltgt b-Pa.
• _ b-Pagt Pb.
• a Pb ltgt a-Pb.

40
Gene Circuit

• a activates a, a activates b, b inhibits a
• First implementation
• a Pa ltgt a-Pa.
• _ a-Pagt Pa.
• b Pa ltgt b-Pa.
• _ b-Pagt Pb.
• a Pb ltgt a-Pb.
• Second implementation exercise to do for next
  week !
• compare both implementations using simulations in
  Biocham

41
Hybrid (Continuous-Discrete) Dynamics

• Gene X activates gene Y but above some threshold
  gene Y inhibits X.
• 0.01X for X gt X Y.
• if Y lt 0.8 then 0.01
• for _ gt X.
• 0.02X for X gt _.
• absent(X).
• absent(Y).

42
Lotka-Voltera Autocatalysis

• 0.3RA for RA gt 2RA.
  0.3RARB for RA RB gt 2RB.
• 0.15RB for RB gt RP. present(RA,0.5).
  present(RB,0.5). absent(RP).