Complex Systems Biology

cf. Life as Complicated System: （current trend)
Enumeration of molecules, processes
detailed models describing the life process
But understanding??
Life as Complex System:
Understand Universal features at a System withmutual dependence between parts and whole
Strategy:1)extension of Dynamical Systems ＆Statistical Physics
Æ Catch consistency between micro-macro levels2) Constructive Approach: (Exp & Theory)` construct simple system to catch universal features'
as simple as possible

Consistency between different levels(1)Cell reproduction vs molecule replication (03-)(2)Reproduction of multicellular organism vs of cells(97(3)Adaptation vs Reproduction(06-)(4)Genetic change vs Phenotypic Fluctuation (03-Constructive Biology Project
theme

experiment

theoryquestion

replicating
systemin vitro
replication withenzymatic reaction

minority
controlorigin ofheredity;
evolvability

replicating cellwith internal
reactions

universal
statistics inreaction dynamics

condition forrecursive
growthcell differentiation.
development

differentiationof E Coilby interaction

emergence ofdifferentiation
rule from dynamics

irreversibility
robustness

Spontaneous
adaptation

Artificial gene network

Adaptive attractor
selection by noise

Robust adaptation
without signalling

evolution

Relevance ofphenotypic fluctuationand dynamics

Genetic assimilationof phenotype
fluct.and dynamics

geno-pheno typerelationshipcellsystemComplex Systems Biology Project (JST,ERATO; KK,Yomo,…)How is recursive production of a cell sustained？
each cell complex reaction networkwith diversity of chemicals;
The number of molecules of each species
not so large

●●●
●●●●
●●●●
●●●
●●●●
●●●●
●●●●

Fluctuations

●●●
●●●●
●●●●
●●●
●●●●
●●●●
●●●●

●●●
●●●●
●●●●
●●
●●●●
●●●●
●●●●

？

Ｔｏｙ Cell Model with Catalytic Reaction Network
ʻCrude but whole cell modelʼ
C.Furusawa & KK、PRL2003
ｋ species of chemicals 、Xo…Xｋ−１
number ---n０ 、n１ … nｋ−１

⚦⢩

model

(Cf. KK&Yomo 94,97):（nutrient）

random catalytic reaction networkwith the path rate pfor the reaction
Ｘi＋Ｘj−＞Ｘk+Xj
some chemicals are penetrable
through the membrane with the
diffusion coefficient D
resource chemicals are thustransformed into impenetrable
chemicals, leading to the growth in
Ｎ＝Σni, when it exceeds Nmax
the cell divides into twocellreaction

:

:
:

:
:
diffusion

catalyze:・・・

K >>1 species

mediumdX1/dt ∝ X0X4; rate equation;
Stochastic model hereIn continuum description, the following rate eqn.,
but we mostly use stochastic simulation

☆Growth speed and fidelity in replication
are maximum at Dc
Dc
1

growth speed
similarity

Growth0.15(only nutrients)
0.6

0.10
0.4
0.05
D = Dc

0.00
0.001

D
No Growth

0.8

0.01
0.1
diffusion coefficient D

0.2

similarity H

growth speed (a.u.)

0.20※similarity is defined from innerproducts of composition vectors
between mother and daughter
cells
・

0
・

Furusawa &KK,2003,PRL
Zipf’s Law is oberved at D = Dc
ni (number of molecules）

rank

number rank
X1 300 5
X2 8000 1
X3 5000 2
X4 700
4
X5 2000 3…….. (for example)

Average number of each chemical ∝ 1/(its rank)-2(distribution of x：ρ（ｘ）∝ｘ ）Confirmed by gene expression data
Human Liver

1e-01

α= -1

1e-02

1e-03

ratio of mRNA

ratio of mRNA

1e-01

1e-04

ratio of mRNA

1e-01

1e+01

1e+02
rank

Human Kidney

1e-02

1e-03

1e-03

1e-04

1e+01

1e+02
rank

1e+03

1e+04

Human Colorectal Cancer
1e-01

1e-02

1e-03

1e-04

1e-04
1e+00

1e-02

1e-05
1e+00

1e+03

ratio of mRNA

1e-05
1e+00

Human Heart

1e+01

1e+02
rank

1e+03

1e+00

1e+01

1e+02
rank

1e+03

1e+04

Mouse ES cell

Mouse Fibroblast Cell
1e-01

1e-02

ratio of mRNA

ratio of mRNA

1e-01

α= -1

1e-03

1e-04

1e-02

1e-03

1e-04
1e-05
1e+00

1e+01

1e+02
rank

1e+03

1e+04

1e+00

C. elegans

1e+01

1e+02
rank

1e+03

Yeast
1e-01

ratio of mRNA

ratio of mRNA

1e-02

1e-03

1e-04

1e-05
1e+00

1e+01

1e+02
rank

1e+03

1e+04

1e-02

1e-03

1e-04
1e+00

Later confirmed by several other groups

1e+01

1e+02
rank

1e+03

Formation of cascade catalytic reaction★ 栄養から直接生成される成分の数を多くなるようパラメータを変えて、階層構造が判りやすくした例
100000ni

Catalyzechemicals ofhigher rank
mainly

４

10000

３

1000

２

100

１

10

Rank of ni
With conservation law,
The exponent -1 is explained
Mean-field type (self-consistent) calc.)

1
1

10

100

1000

10000

100000

１：minority molecules
２：catalyzed by 1, synthesized by resource
３：catalyzed by 2：

• Remarks:（０）Universality(1) Evolution to the critical state (with Zipf law) is
confirmed numerically(2) Evolution to scale-free network appears lateras embedding of power-law abundances into
network (Furusawa,KK, PRE 2006)（3） Self-organization to critical state, if
transport of ‘nutrition chemicals’ is catalyzed bysome chemicals （no need for choice of D)(instead of simple diffusion) (Furusawa,KK,2007)

Evolution of Network to satisfy Zipf’s law?

Yes

Critical D value depends on connectivity in the network;
mutation of network + selection Æ approaches Zipf’s law

Furusawa

Zipf’s law holds, irrespective of network structure, but
Later, the connectivity in the network approaches“scale-free” network through evolution.
statistical properties; embedded into network structure
Dynamics (abundance) first, structure (equation for dynamics) later
evolutionary embeddingof dynamics into network
４
2

probability for a path tochemical with abundances x
is selected; q(x) Ætransformation of abundance
distrb. to connectivity distrib.
Furusawa、KK, PRE, 2006

evolvedinitialSo far average quantity of all components;
Next question: fluctuation by cells:
distribution of each Ni by cells
Log normal distribution !

Furusawa,..
KK,
Biophysics2005
e.g.
cell1 X1 10000
cell2
8000
cell3
15000
cell4
20000
…..
histograrm

Each color
gives
different
chemicalspecies

LOG SCALE

Experiment; protein abundances measured by fluorescence+flow-cytometryLog-normal
Distribution
Confirmed
experimentally

Furusawa,Kashiwagi,,Yomo,KK

Also studied in GFP synthesis in liposome

☆Heuristic explanation of log-normal distribution
Consider the case that a component X is catalyzed byother component A, and replicate; the number --NX、NA
d NX /dt = NX NA
thend log( NX )/dt = NA
If、 NA fluctuates around its mean < NA＞, with fluct. η（ｔ）
d log( NX )/dt = ＜NA＞ ＋ η（ｔ）
log( NX ) shows Brownian motion Æ NX log-normal distribution
too, simplified, since no direct self-replication exists here
But with cascade catalytic reactions, fluctuations are

successively multiplied, (cf addition in central limit
theorem.);Hence after logarithm, central limit th. applied☆Cascade leads to multiplicative propagation ofnoise (at critical region)
Ａ

Ｂ

Ｃ

Ｄ

Ｅ

Propagation of fluctuation, feedback toitself, leading to log-normal distribution tail.
Cf. If parallel,Ａ
Fluctuations come in parallel:
Usual central limit theorem is valid;normal distribution.

Replicating artificial cell (experiment)(ÅÆtheory; fluctuation, minority control)translation

replication

division

RNA polymerase
RNA polymerase geneRNA（Yomo‘s group)

By 菅原ら(東大総合分化）(Sugawara’s group)

Tranlation in liposome
Continouos division of liposomes
RNAreplication in liposome

Questions(1) All chemicals have such large fluctuations?
Important ones are protected??
Origin of heredity (genetic information)
why is there genotype and phenotype
Minority control mechanism (discreteness important)(KK,Yomo JtheorBiol.2002)(2)Large phenotypic fluctuation Æ
relevance to biology ?
ans. evolution (Sato et al., PNAS, 2003) adaptation,….
--- recall in standard evolutionary genetics,
only the distribution of gene is discussed, by assuming
unique phenotype from a given genotype• Phenotypic Fluctuation Æ Relationship to Evolution?
selection is based on phenotype（activity, size, protein abundances, fluorescence,…),butin standard evolutionary genetics;
gene a Æ phenotype x uniquely determinedmostly discusses the phenotype distribution
as a result of genetic variation
Phenotypic fluctuation of isogenic organisms
ÆP(x; a) x—phenotype, a – gene-------------------------------------------------------------------

Artificial selection experiment with bacteria
Selection to increase the fluorescence of protein in bacteria

Schematic drawing of selection process〜2,000 clones
Eyes

Mutagenesis

1st screening

・・・・・ 〜30 clones
Spectrofluorometer
5〜8 clones

2nd screening

Spectrofluorometer
The highest clone

FACS analysis

Ito,Yomo,..

Sato,Ito,
Yomo,KKPNAS(2003)

Fluctuation ---- Variance of phenotype of clone
Organisms with larger phenotypic fluctuation higher evolution
speed;- change of phenotype per generation per mutation -``Response against mutation+selection’’ResponseÅÆ Fluctuation

So-called fluctuation-dissipation theorem in physics:
Force to change a variable x;
response ratio = (shift of x ) / force
fluctuation of x (without force)
response ratio proportional to fluctuation
originated by Einsteinʼs paper a century ago…
Generalization::(mathematical formulation)
response ratio of some variable x against the changeof parameter a versus fluctuation of xP(x;a) x variable, a: control parameterchange of the parameter a Æpeak of P(x;a) ( i.e.,<x>average ) shifts< x > a +∆a − < x > a∝< (δ x) 2 > a =< ( x − < x >) 2 >∆a

Fluctuation-response relationship (generalized form)
Gaussian distribution of x; under the parameter a
at a=a0
Change the parameter from a0 to a

Artificial selection experiment with bacteria
for enzyme with higher catalytic activity
for some protein with higher function
Change in gene (parameter; a) ⇒``Response’’ ------ change of phenotype <x>(e.g.,fluorescence intensity)
per generation per (synonymous) mutation rate
Fluctuation ---- Variance of phenotype x of clone
Fluctuation in the phenotype x of clone⇔ speed of evolution to increase <x>(proportional or correlated)Naïve expectation:
Just propt to mutation rate

Fluctuation-response relation
Phenotype fluct. × mutation rate

(Evolution Speed per generation)
Sato,Ito,Yomo,KK, PNAS 2003

• Confirmation by numerical evolution experimentby the reaction-net cell model
Mutate the network (‘gene’) with mutation rate μ, (rewire the pathof the network with the rate) and select such network
having highest concentration ｃ of a specific chemical
1. Prepare initial mother cells.
2. From each parent cell, mutant cells
are generated by randomly replacing
reaction paths,
with mutation
rate μ3. reaction dynamics of all mutants are
simulated to determine phenotype x
phenotype x = log (ns)

4. Top 5% cells with regard tophenotype x are selected as parent
cells of next generation

Confirmation of Fluctuation Dissipation
Theorem by reaction-network cell modelFluctuation of x=log c

μ=0.01
0.03.0.05

Increase in average x

Furusawa,KK 2005

（１）the use of log(fluorescence), because
log x is close to Gaussian distribution inexperiments(２) New mystery？ phenotype fluctuation ofclone vs evolution speed in contrast to
evolution speed ∝ phenotypic fluctuation by
genetic variation（Vg): (fundamental theorem ofnatural selection; established)pheno fluct of clone Vp∝ pheno fluct by gene variation Vg？
（fluct by noise ∝ variation in ‘equation’)Follow the spirit of Einstein；micro-macro consistencyÆBrownian motion

Consider 2-variable distrb
P(x=phenotype,a=genotype) =exp(-V(x,a))
Keep a single-peak (stability condition).
KK,Furusawa, 2006 JTB

Hessian condition

Up to this point pheno(x) and geno (a) aretreated in the same way.
Then given a, the peak(average) phenotype isx0(a)--function of a --

Gene a is also regarded as variable: Now, we discussstability P(x,a) as 2 variables instead of P(x;a).
write as P(x,a)=exp(-V(x,a))
Condition if it keeps a single-peak solution (stability
problem as in thermodynamics).

Hessian condition

Up to this point pheno (x) and geno (a) are treated in
the same way. Then given a, the peak (average)
phenotype is x0(a)--function of a --x0(a); the average phenotype x of a clone ofgenotype a; x is given as a function of a
Hessian Æ
Vg ≦ Viｐ
Vg= Variance of x, due to the distribution of gene a

Vip＝Variance of x, for the isogenic organism(i.e., a= a0,fixed)
Vg increases with mutation rate μ, thus, there is criticalmutation rate μｃ beyond which the distribution
becomes flat (error catastrophe).
As Vg ∝μ approximately,
Vｇ = (μ／μc) Viｐ
Since Vg ∝ evolution speed （Fisher theorem)
evolution speed ∝ μ Viｐ Æ fluct-response relation

We can do the analysis by using Gaussian 2-body distributionfunction for phenotype x and gene a; around a=a0, and
x=X0;, with coupling between x and a (variance of a isthe mutation rate μ

Stability condition

For high mutation rate single-peakis not sustained :error catastrophe?

Change of distribution
through evolution
Distribution of phenotype
x of a clone
Æ Vp

Log(concentration)

Distribution of phenotype
x over mutants (genetic
variation)
Æ Vg

Phenotype fluct. (Vp) vs Gene Fluct. (Vg) in the evolution of
toy cell model
Vp: fluct. for given network, Vg: fluct. by network variation
μ

Vg

〜μmax

μ

Vｐ=Vg

Vp
Phenotype fluctuation of clone

variance of log(x),x is the concentration of the molecule
Beyond Darwin with the spirit of Einstein!

• ??? to the theory• P(x,a) rather than conditional probability (TRICK)‘’Genetic-Phenotyic correpondence’’what phenotype can vary ÅÆ
what gene can change
fluctuation of variable (micro) vsvariation of equation (genetic evolution)(cf Waddington’s genetic assimilation)
Q: Why error catastrophe when Vg>Vip?
Robust evolution is possible only under noise-counterintuitive ;it says phenotype noise isimportantÆ gene-net modelA simple model for Geno-Pheno relationship;
Model:Gene-net(dynamics of stochastic gene
expression ) Æ on/off state
Xi – expression of gene i : + (on) – (off)Gaussian white
M;total number of genes, ｋ: output genes
Noise strength σ• TaskStarting from -1,-1,-1,,-1(all off)
xi i=1、2、‥・・、k are +1(on) （Target Gene Pattern)
Fitness

F= −（Average number of off x_i）

Genetic Algorithm
Select networks with higher <F>
top--<F>=0
Choose top n networks
among total N,and mutate
with rate μto produce N networks（μ:fixed mutation rate)Result of evolution
Top:reaches the fittest

σfaster for lower noise(σ)

Lowest; cannot evolvefor low noise(σ)

LowestTopσ＜σｃFitness Distribution
σ＜σｃ --low fitness mutants distributed
σ＞σｃ −
eliminatedthrough evolution

Existence of critical noise level σcbelow which low-fitness mutants accumulate(error catastrophe)(1) Vip≧Vg forσ≧σcVg(2)Vg→Vipas σ→σcσcVp=Vg(3) evolution progresses
only for Vip ≧Vg
Vp(4) Vip∝Vgthrough evolution
courseTheory confirmed
KK,PLosOne,2007

Why?; difference in basin structureσ＞σc Æ large basin for target attractor(robust, Δ（distance to basin boudary) ↑σ＜σｃ Æ only tiny basin around target orbit
Δ remains small
Basin Volume for
Each fitnessÆGlobal constraint to potential landscape(funnel?)why threshold?
choose paths to avoid turning
pts within σ (noise)Δ

Mutation→ touches turning
points within range of μΔsmall σ −＞an orbit with small Δ
can reach the target
Δ
Δ

Deviation of basinboundary (turning points)
by Noise −＞δpby Mutation -> δｇ
Vg ~ (δｇ/ Δ）^2
Vip ~（δｐ/Δ)^2
Δ increases
ーー＞robustness
increases
if δｇ＞δｐ,
Δ~distance to turning points
mutation destroys(basin boundary)
the history
ÆVip>Vg necessaryfor evolution of robustness

Discussion• Developmental Robustness to noise ---- Vip• Robustness to mutation in evolution ----Vg
When Vip>Vg, both decrease, i.e., robustness
Noise is necessary for evolution of robustnessVip ∝ Vg ÆDevelopmental robustness and genetic(evolutionary) robustness are linked• Generality of our result; For a system satisfying:(1) fitness is determined after developmental dynamics（2）developmental dynamics is complex(eg.,with distributed catastrophic pt leading to error(3) effective equivalence between mutations and noisewith regards to the consequence to fitness(Æ genetic assimilation by Waddington)

Spontaneous Adaptation• For all possible changes in environment,
signal transduction network is alreadyprovided?
• Or, is there any general (primitive)
mechanism to make spontaneous
adaptation?
• Æ Constructive Experiment with artificial
Gene and theory assuming only growth
condition and stochsticity

(ex) Adaptive response without signal transduction
Unexpected； beyond designed
Embedded gene network
Selection of preferable state
Phenomenological theory of attractor selection
rfp

lacI

tetR

gfp

folA

glnA

Ptet

Env. Without glutamine

Mutual inhibition

Gultamine synthetase
Ptrc

Ptrc
rfp

tetR

folA

lacI

gfp

glnA

Ptrc
rfp

tetR

folA

lacI

gfp

glnA

Ptet
Ptet

Rich environment

Enzyme for Tet

Env. Without Tetrahydro..

fluctuation
Metabolic activity
Kashiwagi,Yomo

Theory of attractor selection by
activity and noise• Embedded network: each of the two can be selected
equally. However, ‘good’ attractor in eachenvironment is selected. Why?
• Due to hidden signal network?
NO!: verified by exchanging the promoter• After each state is attracted with 50%,
cells in a ‘bad’ attractor cannot grow,cells in a good attractor can grow, so that
good attractors are selected?
NO!; the process occurs without (or before) the cell
division process
Novel Mechanism of Spontaneous Adaptation (without
the use of signal transduction) should exist!

• Basic Logic
dx_i/dt=F(Activity）f(x_i)-G(Activity)x_i+η(t)F,G: increase with activity.
active: synthesis, degradation both are fastη Æ noise
Active state ： both Ff and Gx are large
deterministic part >> noise
Poor state ： both Ff and Gx are small
deterministic part ~ noise
Switch from Poor state to Active state by noise(Kashiwagi,Urabe,kk,Yomo; submitted)

The mechanism for adaptive response by attractor selectiondsyn(act )− deg(act ) × m1 + η1
m1 =2
1 + m2dtdsyn(act )− deg(act ) × m2 + η 2
m2 =2
1 + m1dt

6act
syn(act ) =;deg(act ) = act ;2 + act

d
pro
act =
− cons × act
NutthreadNutthread
_
_
dt
1 n1
2
((
) + 1) × (() n2 + 1)
m1 + Nutrient1
m2 + Nutrient 2

m1

m2

Adaptive Response of thegenetic network to aenvironmental changeact

No Nutrient1

No Nutrient
Nutrient2
No depletion

No depletiondsyn(act )− deg(act ) × m1 + η1
m1 =
2
dt
1 + m2dsyn(act )− deg(act ) × m2 + η2
m2 =
2
dt
1 + m16act
syn(act ) =;deg(act ) = act ;2 + act

1.2

AttractorW

1

act

0.8

0.6

Attractor1

0.4

d
pro
act =
−cons×act
NutthreadNutthread
_
_
dt
1 n1
2 n2
((
) +1)×(() +1)
m1+Nutrient1
m2+Nutrient2

0.2

0
0

0.5

1

1.5

2

2.5

m1

Fraction that adaptive attractor is selected3dm1 = 0dt
2.5

Nutrient1 = 0
Nutrient 2 = 10

m2

2

1.5

separatrix

1dm2 = 0dt

Noise strength
0.5

Adaptive Attractor
0.5

1

1.5

m1

2

2.5

3

3

• Growth-Induced-Attractor-Selection (Furusawa kk)• Basic Logic
dx_i/dt=f(x_i)-S({x_j})x_i+η(t)SÆ dilution effectη Æ noise
Active state ： both f and S are large
deterministic part >> noise
Poor state ： both f and S are small
deterministic part ~ noise
Switch from Poor state to Active state by noise
Selection before reproduction
General logic in a system with growth andfluctuation

Gene network -> a huge number of attractors coexist with different
growth speeds

Spontaneous selection of optimal growth states
General in a system with noise and growthSummary
Consistency Principle for Biology-- replication of molecules and cells :Universal Laws(-- replication of cells and cell ensembles)--- genetic and phenotypic changes(+speciation)--- adaptation of internal cellular state and growth・Biological relevance of phenotype fluctuations?
ÆPhenotypic Fluctuation ∝ Evolution Speed
Æ Relation between(isogenic)phenotype fluctuationｖｓ
phenotype variation by mutation• Robustness to mutation and to developmental noise
are linked• Growth system Æ general adaptation by noise• consequence of steady growth system

Collaborators
Chikara Furusawa
Katsuhiko Satoexperiment
Tetsuya Yomo
Yochiro Ito
Akiko Kashiwagi
Most papers (biology,
Dynamical systems)
Available at
http://chaos.c.u-tokyo.ac.jp
ERATO Complex Systems Biology Project(2006,August)