######### ANALYSIS AND SIMULATION OF CELL NETWORKS

+-------------------------------------------+--------------------------+
| **Principal Investigator**                | Phone: 630-252-5941      |
|                                           |                          |
| Evgeni Selkov, Senior Computational       | Fax: 630-252-5986        |
| Biologist                                 |                          |
|                                           | Email:                   |
| Mathematics and Computer Science Division | <selkovsr@mcs.anl.gov>   |
|                                           |                          |
| Argonne National Laboratory               |                          |
|                                           |                          |
| 9700 South Cass Avenue                    |                          |
|                                           |                          |
| Argonne, IL 60439                         |                          |
+-------------------------------------------+--------------------------+

### 

### The goal of the proposed effort in analysis and simulation of cell networks is to increase our understanding of life phenomena at the level of bacterial cells. To achieve this goal, we propose a two-pronged effort focusing on (1) development and analysis of large-scale mathematical models of selected sequenced photo- and heterotrophic bacteria and their populations and (2) fitting of the models to quantitative experimental data.

### Specifically, we plan the following developments in mathematical modeling:

### We will automate development, debugging, solving, optimizing, analyzing, and visualizing mathematical models of large cell networks.  The complexity of models will range from stoichiometric matrices to nonlinear ODE and PDE equation sets with possible stochastic perturbations all covering large parts of cells or entire cells. 

### From the available sequenced genomes, we will compute a universal metabolic core, common for all sequenced autotrophic bacteria, and use it as a basis for the initial mathematical simulation of the central metabolism, common for such bacteria.  

### We will extend this core flux model to a model with specific features of *Synechococcus* WH8102, SFM, common for other cyanobacteria, and to a model of *Pseudomonas fluorescens*, PFM, common for other organotrophs.  

### Based on proteomic data from *Pseudomonas fluorescens* and *Synechococcus*, we will analyze the metabolic differences between plaktonic and adhered (biofilm) bacteria to develop a more complete characterization of the differences in bacterial behavior between these two environments.

### ***Background and Significance***

It is well known that an open biochemical reaction system described by a
deterministic model of three nonlinear ordinary differential equations,
ODEs, can generate a chaotic behavior \[Goldbetter 1996, Gurel, 1983,
Sel'kov 1980\]. There is no way to precisely describe or explain such
complex behavior by using a spoken language or quite scientific
biological terminology. Language of mathematics and mathematical models
are the only way toward precise description and deep understanding of
complex biological systems.

Therefore, once a genome is sequenced and the encoded metabolism is
reconstructed, building a mathematical model of the metabolism to
understand how it works is the next obvious step.

Serious obstacles, however, make simulation of the metabolic networks a
challenging task. The most significant of these are the following:

1.  *Lack of information on regulatory mechanisms*. In spite of a very
    successful use of the metabolic reconstruction technology in
    academia and in industry, the current applications of it are limited
    by predicting a metabolic skeleton of the cellular organization and
    chromosomal clustering of functionally related genes. The whole
    complex network of regulatory interactions is currently beyond
    resolution of this technology.

2.  *Lack of numeric data*. For a realistic quantitative dynamic model
    of the whole cellular metabolism, one needs an enormous amount of
    numeric data on enzyme and transporters, their activities and
    substrate specificities, kinetic and equilibrium constants,
    concentrations of proteins and intermediates, and so forth. Even for
    a minimum model of the smallest reconstructed genome of *Mycoplasma
    genitalium*, one may require hundreds of such numbers. It would take
    many years to get the required data experimentally.

3.  *Insufficiency of data for a single cell*. Simulation of the dynamic
    behavior of a single cell progressing along its cell cycle requires
    qualitative and quantitative data on the temporal behavior of many
    components of a single cell. Measuring a dozen of parameters along
    several cycles of a synchronously growing bacterial culture is a
    heroic endeavor that is no longer attempted. It is not well known
    that intercellular interactions in an apparently synchronous and
    homogenous cell culture can lead to a wide spectrum of dynamic
    behavior that can be attributed to a single statistically averaged
    cell*.*

4.  *Lack of knowledge about the functions of many proteins*. It is
    common to point out that the essential percentage of every sequenced
    genome constitutes hypothetical proteins. Unfortunately, there is no
    indication that all of them are of even minor importance for the
    cellular life.

5.  *Uncertainty about intracellular conditions*. Intracellular
    conditions may be quite different from those used in experiments in
    vitro to determine the kinetics of purified enzymes. The debate has
    continued for decades.

6.  *Complexity of the mathematical model.* The mathematical model may
    be too complex to be a good substitute for a real cell. Even the
    simplest mathematical dynamic model of growing *M. genitalium* may
    have at least tens of variables and many tens of parameters. A
    complete bifurcation analysis of such model may seem to be barely
    possible. Manual encoding and debugging mathematical models derived
    from genomes much larger that the *M. genita*lium represent a tough
    technical problem even if all other problems are solved.

7.  *Nomenclature problems*. All parts of the living cell (mini-, midi-,
    and macromolecules, their aggregates, functional blocks, control
    signals, mechanisms, locations, and many others) have names.
    Unfortunately, there is no universal nomenclature─with the rare
    exception of chemical nomenclature and nomenclature of enzymes─and
    even these two nomenclatures are only recommendations that are
    ignored by many scientists.

Mathematics of complex biological system developed during the past
several decades has found many techniques to address these problems. We
will use these techniques in our mathematical simulation of
reconstructed cellular networks.

##### Preliminary Results

**Occam's razor and other theoretical methods of reducing complexity**.
Facing an enormous complexity of biological organization, theoreticians
must find the ways of making complex system much simpler and
transparent. The well-known economy principle of William Occam
\[Thorburn, 1915\] remains among the most fundamental philosophical
ideas underlying the formulation of the most successful model of
biological systems.

Hodgkin and Huxley \[Hodgkin, 1952\] neglected all details of excitable
plasma membrane of giant axon of *Loligo* (molecular structure of its
ion channels, slow regulation of channels and interaction with the axon
metabolism, etc.) and developed an elegant quantitative mathematical
model of membrane excitation. Their model that not only explained the
observed phenomena but predicted a number of generic channel features
that were confirmed decades later when the 3D structure of
voltage-controlled ion channels became known.

Mathematics has developed powerful tools of systematic dramatic reducing
complexity of mathematical models. Many of them basically exploit a very
simple idea of neglecting small numbers when they are compared with
large ones. Among these mathematical tools is Tikhonov's theorem
\[Tikhonov, 1952\] and its generalizations \[Vasilieva, 1967\]. The
general idea of this theorem is that if a mathematical model of a
dynamic system can be represented as an ODE set

***x*'**= ***f(x, y),***

***y'= g(x, y)/**e**,***

where *x* and *y* are vectors of slow ***x*** and fast ***y***
variables, ***f*** and ***g*** are vector functions, and *e* \<\< 1 a
small scalar number (a famous "small parameter"), then the conditions
specified by the theorem allow a limit transition e  0 that gives a
much simpler asymptotically perturbed system

***x*'**= ***f(x, y),***

***g(x, y) = 0,***

solutions of which are very close to the solutions of the original one.
The beauty of this reduction is that the details of the fast motion can
be treated separately in the short time t~s~ = e\*t, t is the original
time. Rewriting the original model for the short time t~s~ gives it a
new look

***x*'**= e***f(x, y),***

***y'= g(x, y).***

The model is now to be used to analyze details of the fast motion. This
form suggests making another radical simplification of the original
model by a new limit transition e  0. This transition friezes the slow
variables (***x*'**=0, ***x***= *const*) and gives a simplified version
of model

***y'= g(x, y),***

describing fast converging of variables ***y*** to their
quasi-stationary values determined by the solutions of ***g(x, y) =***
0.

Tikhonov's theorem is probably the most frequently used mathematical
tool used in simulation of complex cellular networks (see general
references in \[Reich, 1975, 1981\]). It is suffice to say, that all
rate laws of enzyme reactions, published in the literature since the
classical article of Michaelis and Menten \[Michaelis, 1913\] are
derived as quasi-stationary approximations justified by the Tikhonov's
theorem.

**Mathematical analysis and simulation of cell energy metabolism.** In a
series of publications, Selkov and his coworkers showed \[Selkov, 1979;
Reich, 1981\] that an uncontrolled stoichiometric backbone of
intermediary metabolism is capable of fine stabilizing of relative ATP
concentration (i.e., Atkinson's adenylate charge, or *energy charge*
\[Atkinson, 1969\]) against variation of metabolic demand in a wide
range \[Ivanitsky, 1978; Selkov, 1975, 1979; Reich, 1981\]. This
theoretical prediction was later confirmed in experiments with a
glycolytic system of a cell-free yeast extract treated by an excess of
ammonium ions to desensitize key glycolytic enzymes to any allosteric
signals \[Boiteux, 1984\].

**Two recipies for predicting regulatory mechanisms and missing
paramer.** It has also been shown that stoichiometric structures of
cellular metabolic networks can generate a wide range nonlinear
phenomena (hysteresis, multiple steady states, fine stabilization of
concentrations, self-oscillations, etc.) earlier ascribed to much more
complicated allosterically or genetically controlled networks
\[Ivanitsky, 1978; Selkov, 1975, 1979; Reich, 1981\].

These observations led Selkov \[Ivanitsky, 1978, Selkov, 1979, 1980c4\]
to the idea that evolutionary younger nonstoichiometric regulatory
mechanisms (allosteric and genetic) have evolved during the biology
evolution to improve functioning of stoichiometric networks by
compensating for their functional drawbacks. From this follows the
recipe of theoretical predicting regulatory mechanisms:

-   Formulate a mathematical model of a stoichiometric network.

-   Analyze the model to find the defects.

-   Compute the best control mechanisms and their optimum parameter
    values to compensate the defects.

Selkov \[1980c\] applied this method for predicting possible allosteric
mechanism in purine metabolism. The function of the searched control
mechanisms was to keep absolute concentration of ATP constant. The
literature search confirmed that the predicted mechanisms were
experimentally found in different organisms and tissues, although in
different combinations with each other.

The proposed method helps to solve two heavy problems simultaneously: it
predicts control mechanisms and missing parameter values. It appears
that the method requires one to know the functions of cellular systems
in order to formulate an optimization algorithm. One can argue that such
functions may be completely unknown or our assumptions about them wrong.

There exists, however, a simple way to circumvent this problem. It
requires experimental data on concentrations of metabolites, enzyme
activities, time behavior, and so forth─the more, the better. The data
can be used to search for the best combination of regulatory mechanisms
and parameter values of a system model satisfying the best fit to the
data. This variant of optimization does not require any information a
priory about system functions.

In fact, we do assume one formal function -- keeping the model behavior
as close as possible to the experimentally observed one. Yet, we are
free to formulate our hypothesis about the possible functions to
optimize controls and parameter them. We thus have a regular way of
testing our hypothesis about possible system functions.

We will use such optimizations to extend the limits of the
reconstruction technology and to build comprehensive qualitative models
of bacterial cells. The MCS Division has considerable expertise and
excellent software, and hardware to support this part of our project

**Studies of biochemical oscillations and the cell clock**. Studies of
biochemical oscillations in metabolic and genetically controlled systems
were very popular in the 60s and 70s; see references in \[Goldbetter,
1996; Gurel, 1983\]. These studies were always heated by a hope to
understand mechanisms of the cell clock \[Bunning, 1972; Edmunds, 1984,
1988\]. The cell clock is a hypothetical self-oscillatory mechanism
controlling periodic changes in the cell leading to cell divisions.

Selkov \[Sel'kov, 1970\] has suggested that the cell clock is a limit
cycle phenomenon related to thiol metabolism. His very speculative model
assumed that the leading factor generating auto-oscillations in thiol
content is a positive feedback in formation of reduced thiols. The model
has been developed to theoretically explain experimentally observed
oscillations in thiol content (aka Rapkine's cycle \[Rapkine, 1931\]) in
many experimentally studied organisms and in cell-free extracts of sea
urchin eggs \[Mano, 1968, 1969, 1970\]. Later, Mano extended his studies
on thiol clock and has shown that the thiol oscillations propagate
through the cellular metabolism of cleaving eggs and cause cyclic
response in protein synthesis, RNA polymerase, and DNA polymerase
\[Mano, 1971-1977, Mano, 1976\]. Thus, Mano confirmed Selkov's idea
\[Selkov, 1970\] that the thiol clock controls the entire cell
metabolism.

Tin spite of these successes, the purely qualitative model of Selkov
\[Sel'kov, 1970\] could not give answers to the following basic
questions::

1.  What is the primary function of the clock if the clock ticks with a
    period much shorter or longer than a day?

2.  How can the clock generate very slow oscillations with the period of
    about a day in the metabolism that have characteristic times of
    several minutes, tens of minutes, or even hours?

3.  What does make the Circadian version of the clock so stable?

4.  What are relationships between the clock and the cell division

5.  What is the precise molecular mechanism of the clock?

6.  Is it universal, or there are many versions of the clock?

Analysis of stoichiometric models of cell energy metabolism suddenly led
Selkov and coworkers \[Ivanitsky, 1978, Selkov, 1979, 1980c; Avseenko,
1980, Boiteux, 1980\] to a discovery of a critical role of temporal
separation of futile cycles of cell energy metabolism. The model
computations show that an uncontrolled substrate cycling in the futile
cycle formed by two cytosolic enzymes 6-phosphofructokinase, PFK, (EC
2.7.1.11) and fructose-bisphosphatase, FBPase, (EC 3.1.3.11) leads to
the negative net production of ATP. There is no way to prevent this
wasteful substrate cycling besides a temporal organization.

Such organization requires an auto-oscillator (a clock) to control these
two enzyme activities reciprocally. The theory predicted that the
oscillator must have the oscillation period as long as possible to
minimize the futile recycling during the transitions between two
alternating states \[Selkov, 1979a, 1979c\]. Minimization of futile
recycling requires the transition times as short as possible. From this
point of view sine waveform is among the worst possible and must be
avoided \[Selkov, 1979c\]. Such oscillations with a very short
transition time can be generated by a mechanism of product activation
only.

This conclusion about the required waveform kills all theories about any
cell clock mechanism based on a negative feedback \[Goldbetter, 1966\].
The negative feedback oscillators require at least three consequtive
reactions \[Selkov 1967, Morales, 1967, Aponin, 1971, Hunding 1974,
Othmer, 1976, Heirich, 1996\]. These three steps act as a low band pass
filter and cutoff effectively all high frequencies except the frequency
of excitation for which the phase delay is 180 degrees. This makes the
waveform, generated by negative feedback oscillators, quasi-sinusoidal.
The smooth oscillations observed so far in many experiments with cell
clocks of different organisms can be explained by a wide phase
distribution of individual cell clocks in a not ideally synchronous cell
population. It is clear that one needs a special research to study
population effects on the waveform generated in many individual cells.
Without such study, any effort to guess about a possible waveform of
intracellular oscillations will have a very low chance to reflect
reality.

Selkov \[Selkov, 1972, 1979, 1980a; Schulmeister, 1978\] found the
answer to the problem of the long oscillation period. He determined that
a fast reversible exchange between a storage compound like glycogen with
the glycolytic oscillator makes the oscillation period very long. His
estimates made specifically for a glycolytic oscillator \[Selkov,
1980a\] gave oscillation period about a day.

The prediction that glycogen must be a participant of the Circadian
clock was confirmed in experiments with *Gonyaulax polyedra* \[Dunlap,
1980; Rensing, 1980\]. The work showed that phase shifts of the clock
induced by cycloheximide are connected with the corresponding phase
shits in the glycogen content.

The authors explain the phase shits by an inhibitory action of
cycloheximide on protein synthesis. However, the concentration of the
inhibitor was very low and inhibited protein synthesis only partially.
The hypothesis that the inhibitor acts specifically on synthesis of some
clock protein(s) is not convincing.

In rats, cycloheximide activates a quick disappearance of liver glycogen
which cannot be explained be action of the inhibitor on protein
synthesis. If the inhibitor had the similar effect on *G. polyedra* the
addition of it in the phase of glycogen accumulation would result in
phase delay. Addition of the inhibitor in the phase of glycogen
degradation should result in the phase advancement. The authors
\[Dunlap, 1980; Rensing, 1980\] described precisely that kind of action
of cycloheximide.

Glycogen plays a central role in Circadian temporal organization of
cyanobacteria. Direct measurements of glycogen content in synchronous
culture demonstrate large a Circadian oscillations in glycogen content
with the expected wave form \[Schneegurt, 1994\]. Thus, participation of
glycogen in the Circadian clock mechanism seems well documented.

Selkov \[1979a\] predicted that multiple negative feedbacks protecting
all aspects of cellular metabolism should stabilize the clock period.
This prediction has been tested with a large mathematical model
describing the cell energy metabolism \[Avseenko, 1987; Selkov, 1989\].
The model included glycogen metabolism connected with
oscillation-generating futile cycle, catalyzed by 6-phosphofructokinase
(EC 2.7.1.11) and fructose-bisphosphatase (EC 3.1.3.11), and four well
known negative feedbacks controlling glycolysis/gluconeogenesis.

This model demonstrated all main features theoretically expected from
the Circadian clock. The model generated very stabile Circadian
oscillations. The oscillations effectively suppressed wasteful substrate
cycling in three major futile cycles of the metabolism. The temporally
organized energy metabolism maintained a high and stabile level of ATP
to power an external metabolic load. The optimized model had a very
large valley in its parametric space -- in this valley, variation of
parameters in a wide range did not resulted in a substantial variation
of oscillation period. Thus, it has been proven that the temporal
organization of energy metabolism can behave as a Circadian clock with
all its postulated attributes.

This model, however, did not solve problem of the molecular mechanism
yet. The oscillation mechanism of the model had two positive feedbacks,
activation of PFK by FBP and reciprocal inhibition of FBPase by FBP.
These positive feedbacks are missing in cyanobacteria. Therefore, we
need to find functional analogs of these positive feedbacks.

Our preliminary study, based on the comparative analysis of sequenced
and reconstructed genomes, predicts that the most probable positive
feedback mechanism will include reciprocal control of several key
enzymes by a thioredoxins, the redox status of which will be tightly
controlled by an intermediate of the lower part of the glycolytic
system. Such a model could display a strong coupling between energy
metabolism and macromolecular metabolism. It will show phase shits to
induced by addition of redox dyes, intermediates like pyruvate, lactate,
inhibitors of protein synthesis, and respiratory inhibitors.

The problem 5 listed above has been more or less solved experimentally.
It is known that in cyanobacteria the clock ticks in stationary cells.
This means the clock is independent from division cycle but gates the
division.

**Enzyme reactions of real cells are open twice**. They are open with
respect to their reactants and products and with respect to enzyme
molecules as well. A slow spontaneous or catalytic degradation of enzyme
molecules in is usually compensated by a *de novo* biosynthesis to keep
a certain steady state level of enzymes. Selkov and Nazarenko \[Sel'kov
1980, Nazarenko 1983\] have shown that a double open enzyme reaction
behaves quite unexpectedly \-- it can generate oscillations with a very
long oscillation period T \~ 1/k~E~ where k~E~ is enzyme degradation
rate constant. In a simple single-reaction model oscillations are slowly
converge to a stable steady state. There are many ways of converting
such oscillations into limit-cycle oscillations.

In particular, Brownian noise due to a low concentration of enzymes may
continuously perturb the slowly damping oscillations. This eventually
will produce a limit cycle the size of which will be strongly dependent
on the average amplitude of the noise. Such transformation of a stable
focus into a limit cycle has been studied in with simple variant of a
glycolytic model \[Ebeling, 1983, 1985\]

**Derivation of rate lows for complex controlled reaction mechanisms.**
The founders of modern theories on mechanism of action of allosteric
enzymes \[Monod, 1963, Koshland, 1966, Cornish-Bowden ,1995\] did not
bother themselves by the complexity of real enzyme reaction. Their
models describe the simplest possible enzyme reaction catalyzed by
active sites of oligomeric enzymes, that is an irreversible
transformation of a single substrate into a product.

The vast majority of enzyme reactions are reversible and they have
two-three substrates and two-three products \[Enzyme Nomenclature,
1992\]. Multiple attempts to derive rate laws for such reaction
catalyzed by oligomeric enzymes failed before Popova and Selkov
\[Popova, 1975, 1976, 1978, 1979; Ivanitsky, 1978\] have found a regular
way of deriving such complex rate laws and demonstrated that a
generalized multisubstrate model gives a good fit to
6-phosphofructokinase and pyruvate kinase reactions from E. coli
\[Malkova, 1980, Miakishev, 1983\]

### **Experiment Design and Methods** 

Our simulation strategy will be based on gradually increasing the
complexity of mathematical models and the levels of biological hierarchy
matched by the simulation. It has been shown \[Edwards, 2000; Cowert,
2001\] that such models can successfully represent the metabolic
phenotype of bacterial cells. The flux models later will be upgraded to
kinetic and dynamic models of a single cell. Finally, a series of models
of increasing complexity will be developed to describe the dynamics of
cell populations under different environmental conditions.

Specifically, the steps of developing of steady-state models of
*Synechococcus* and *P. fluorescens* will be the following:

1.  Metabolic reconstruction for available sequenced genomes of
    cyanobacteria and gram-negative organotrophs, phylogenically close
    to *Pseudomonas fluorescens* will be performed. This will require
    analysis of 20-30 phylogenically related genomes for each organism
    to predict conserved chromosomal gene clusters, and therefore to
    increase precision of gene function annotation; use and development
    of standard data type definitions, ontologies, and nomenclatures;
    and updating of the metabolic pathway database. Some assertions of
    functions and pathways will require experimental tests.

2.  Compilation and debugging of the flux core models for
    *Synechocystis* and *P. fluorescens*, based on the metabolic
    reconstructions of these organisms from sequenced genomes.

3.  Development of consistent nomenclature of metabolites, protein
    functions, pathways, compartments, and so forth. Normalization of
    the nomenclatures is a bottleneck in the entire reconstruction
    technology preventing an automation conversion of reconstructions
    into mathematical models.

4.  Diurnal modeling of *Synechococcus* by fitting the flux model *of
    Synechococcus* to two incompatible alternative metabolic states: (1)
    of CO2 fixation and (2) N2 reduction to obtain models of flux
    distribution under these conditions. This step requires experimental
    measurment of metabolic fluxes toward main storage compounds
    (glycogen, cyanophycin, sucrose, trehalose, glycosylglycerol,
    reserve lipids, polyphosphates, etc.) in addition to the main
    metabolic fluxes across the plasma membrane. Experimental measurment
    of the intermediate compounds will be done by using data from a
    chemostat culture culture synchronized by periodic illumination LD
    12:12 hrs.

5.  Fitting the flux model of *Pseudomonas fluorescens*, PFM, to
    experimental data obtained at different modes of batch growth
    (log-phase/stationary phase) and physiological state
    (planktonic/attached)

6.  Sensitivity analysis of these models to allow prediction of the most
    probable control steps in the metabolism affected by different
    physiology states. This analysis will be supported by the cluster
    analysis of gene expression data.

7.  Conversion of the flux models into the simplest stationary kinetic
    model in which all reactions rate laws have no saturation terms.
    This "chemical approximation" will be represented by a large set of
    nonlinear algebraic equations. It will be used for an initial
    evaluation of a large number of rate constant from computed fluxes.
    Optimization software developed by MCS will be used for this
    parameter fitting.

8.  Qualitative and then quantitative analysis of the fitted algebraic
    models to reveal inconsistencies due to lack of regulation component
    in developed "unregulated" model of central metabolism. Such
    inconsistencies will be used in the next step for computing the best
    combination of regulatory mechanisms that leads to a correction of
    the defects. Experimental validation and support will be integral to
    this process.

9.  Formulation of an optimization task to compute the best set of
    allosteric mechanisms and regulatory parameter values that minimize
    the revealed defects. This is a computationally demanding task.
    Initially it will be simplified by a number of technical methods,
    including the power law approximation of control functions and the
    use of blockwise optimization as opposed to global optimization. It
    will be further simplified by assuming that only cross-road
    metabolites are allosteric regulators and that only irreversible
    reactions are potential regulation targets. Experimentally confirmed
    regulatory interactions will be a subject of further detailed
    experimental kinetic studies with purified enzymes.

10. These optimized models will then be converted to more realistic ones
    by introducing saturation terms into rate laws of controlled
    reactions. This will require using detailed kinetic models for
    description of the mechanisms of these reactions. This extension
    will introduce many unknown kinetic parameters, which will be
    computed by fitting the models to experimental data on flux
    distributions.

At this stage the models can be upgraded to dynamic models describing
temporal changes in the cellular metabolism.

#### We propose to address the following six specific issues: (1) building of dynamic models of a single cell, (2)development of the cell clock theory, (3) simulation of synchronous modes, (4) simulation of inhomogeneous cell populations, (5) simulation of possible dynamic behavior in a homogeneous film, and (6) analysis of possible effects of small numbers of molecules on cell population. Below, we discuss each of these in turn.

#### Building of dynamic models of a single cell. This project is focused on the process of adhesion in which bacterial cells, initially free-floating, become adhered and then mature into a structured community called a biofilm. During this process, the physiology of the cell changes dramatically, as does its relationship to its environment. To build a dynamic model of this process requires far more experimental information about the process than we currently have in hand. However, initial modeling will be based on the detailed evaluation of the two states: planktonic (a representative free-floating state) and adhered (a representative biofilm state). The cellular subsystems that need to be shut down and the others that need to be turned on to make this transition will be detailed through a combination of proteomics, metabolic modeling, and analysis of the transcriptome, supplemented by additional experimental data as deemed needed for the detailed modeling. The regulatory pathways will be mapped through a close interaction between our efforts to build a database of signal transduction proteins (Sentra) and our protein interaction mapping work.

#### Development of the cell clock theories and controversies. Because of the different physiologies of cells in planktonic and biofilms, the interaction of the process of adhesion/biofilm formation and the cell clock will be critical for developing a comprehensive understanding of biofilms. This is particularly true in the case of *Synechococcus,* which maintains a cell division cycle of approximately 24 hours that is critical for its incorporation of incompatible metabolic pathways within a single cell.

The controversy about the mechanism of the cell clock involves the
nature of the feedback mechanism controlling it. The popular negative
feedback model is not convincing, however, for a number of reasons.

#### The primary function of the cell clock is to separate in time incompatible biochemical processes.  In all oxygenic N~2~-fixing cyanobacteria C O~2~ fixation coupled with O~2~ evolution is incompatible with N2 reduction. In H~2~-generating *Synechococcus* strains C O~2~2 fixation/O~2~ evolution is incompatible with H~2~ generation by cytosolic hydrogenase.  In both cases, light-driven CO~2~ fixation/O~2~ evolution leads to accumulation of glycogen that is used then to generate reducing equivalents and ATP to power N~2~ fixation by nitrogenase or H~2~ formation by  reversible cytosolic hydrogenase. 

#### Any bioengineering of cyanobacteria to improve their CO~2~ sequestration, nitrogen fixation, H~2~ generation, and biomass production requires optimizing periodic solutions of dynamic models describing the whole metabolism and its regulation.  In nonphototrophic bacteria, the biochemical incompatibility can be not that obvious.  Yet, theoretical analysis (Selkov, 1979; Bioteux, 1980) has shown that any cell energy metabolism containing futile cycles cannot generate net ATP without temporal separation of oppositely directed reactions of such cycles.  In phototrophic organisms, including photrophic bacteria, the oscillation period of the clock is close to 24 hours and very stable.  Such a stable Circadian clock has an additional function of measuring diurnal time.  

#### The temporal organization of incompatible biochemistry requires very fast transitions from one metabolic state to another to prevent overlapping of the incompatible processes (Selkov, 1979).  Rectangular oscillations, typical for positive feedback mechanisms, are the best to meet this requirement.  By contrast, negative feedback mechanisms with their smooth sine wave-like oscillations do not fit the criterion. 

We believe that building the whole cell dynamic model of
*Synechococcus*, a classical organism for studying Circadian cell clock,
and fitting it to experimental data obtained with synchronous cell
populations will be the shortest way of resolving the long theoretical
discussion on the cell clock mechanism.

We will investigate whether the core of the clock mechanism is common
for such different organisms as *Synechococcus* and *Pseudomonas
fluorescens*.  This effort will require building quantitative whole cell
dynamic models of each of these organisms and their comparative
theoretical and experimental analysis.   It will also require extending
the single cell simulation study toward simulation of synchronous
populations of these organisms.

**Simulation of synchronous modes in homogeneous populations of cells
(interacting via common factors of the growth medium.).** Interactions
among cell clocks of individual cells in cell populations can lead to a
synchronization of metabolic and physiology events of individual cells. 
The interactions can be established via common growth substrates,
secreted products, dissolved O~2~ and CO~2~, and specific hormone-like
compounds and toxins. Mutual synchronization of intracellular events may
be vitally important for many species in their natural habitats. 
Obtaining the maximum cell population synchrony is critically important
to enable measurement of intracellular events and tracing their temporal
order.** **

#### We propose to analyze an ODE model of 1,000 *SynechococcusWH8102* cells interacting with each other via different factors common of its cultivation medium including quorum sensing compounds.  We will initially analyze this model under assumption of a homogeneous distribution of cells and ideal stirring of the medium. We expect to carry out a bifurcation analysis to determine parameter domains of existence of different modes of full or partial synchrony under conditions of a complete identity of individual cells and preset parameter scatter. We will determine the conditions favoring a stable auto-synchrony in the population even under essential scatter of individual cell parameters.  We will test hypotheses concerning possible instrumental control mechanisms stabilizing the population synchrony. We will use experience obtained to repeat the analysis with the model developed for *P.  fluorescens.*  

#### Simulation of spatially inhomogeneous cell population. We will convert the ODE models of *Synechococcus* and *P. fluorescens* into a hybrid model of partial differential equations (PDEs), describing diffusion of extracellular compounds into and out of cells, and systems of ODE, each describing a single cell.  Such models will represent a thin 2D layer of 100-100 cells equidistantly or randomly spaced in a square space of cultivation medium. 

#### Similar but much simpler systems have been analyzed, theoretically and experimentally, to study wave propagation in homogeneous chemical and biological media. Extrapolating the results, we expect to find the whole variety of dynamic behavior in such 2D cell communities, including propagation of slow waves of local metabolic synchrony and strong synchronizing effect of periodic illumination for *Synechococcus*. 

#### A more complex variant of the 2D layer of *Synechococcus* cells will also be analyzed.  The layer will model a vertical cross-section of seawater with a gravity force applied along one axis and a free float of the cells.  The buoyant density of the cells will be computed from the current amount of storage compounds (glycogen, cyanophycin, triglycerides, polyphosphates).  The model will describe periodic motion of the cells during the diurnal cycle of illumination.  The maximum content of glycogen by the end of light phase will result in slow sinking the cells.  Possible effects of quorum sensing on this cyclic migration are difficult to predict now. 

#### Similarly, a vertical 2D layer of floating *P. fluorescens* will be analyzed.  Here, the gravity force will act on the cell of the buoyant density dependent on the current content of polyphosphates, glycogen, and lipids.  We expect that the model will predict a number of critical phenomena associated with gradients of O~2~, cells, and temporal changes in the buoyant density of cells: beyond some critical local concentration of O~2~, the cells switch to the fermentation mode associated with accumulating polyphosphates. The increasing buoyant density results in sedimentation of the cells to the bottom.   Simulation of dynamic behavior in a homogeneous film of *Pseudomonas fluorescens.* We will extend the 2D model for *P. fluorescens* by converting it into a 3D model.  The model will describe a cubic cultivation volume with the cells attached to the bottom of the volume as a regular or irregular mesh of cells.  The cultivation fluid can be still or flow with a constant flow along the cell layer. The geometry of biofilm, cell layering and spacing, and so forth will be specified depending on the results of previous simulations and comparisons with experimental studies of biofilms.   Analysis of possible effects of small number of molecules on cell and cell population dynamics. Several investigators have pointed out that the limited number of molecules present in a cell may give rise to stochastic behaviors that must be incorporated in any comprehensive effort at cell modeling. *Synechococcus* and *Pseudomonas* have cytosol volumes of about 1e-15 to 1e-14 l.  With concentrations of metabolites it the range of 0.01 to10 mM, the cytosol would have 1e3 to 1e7 molecules of each molecular species.   A very low concentration of some catalytic or regulatory proteins might be a problem if they interact with each other.  There is no problem, however, if weakly expressing proteins interact with metabolites in the above-mentioned range.  Yet, typical for cytosol compartment pH = 7.0 accounts for only 10 to 100 protons in the compartment.  The ubiquitous presence of protons as a reactant or an effector appears to make the whole metabolism quite stochastic.  

#### These estimates do not take into account the strong buffer effect of the whole network.  Many low-concentration reactants are buffered by their metabolic or binding partners.  For instance, for a single virtual free proton in a compartment at pH = 8.0 there will be excessive buffer capacity of all moieties R, which are in a very quick equilibrium with free protons   H^+^ + R == RH.  With total R+RH of about 10 to 100 mM one could apparently neglect possible very fast fluctuations in a local proton concentration.  We realize, however, that only a detailed computer simulation of a metabolic network running in a cell volume with realistically chosen numbers of molecular species can give a convincing picture about possible stochastic component in the dynamics of the cell metabolism.  We plan to use a dynamic model describing the core metabolism of a phototrophic bacterium to test contribution.

**4. Identification of Classes of Proteins Involved in Biofilm
Formation**

We will identify the following classes of proteins involved in a process
of biofilm formation in *Synechococcus* WH8102 and *P. fluorescens*:

-   Proteins involved in quorum sensing (lux-genes, signal transduction
    proteins) using the Sentra database technology based on domain and
    motif analysis

-   LPS biosynthetic genes in *P. aeroginosa* and *P. fluorescens* (upon
    availability of a genome) will be predicted using gene function
    predictions algorithms and tools described in the Bioinformatics
    section of this proposal

-   Secreted proteins (e.g., proteases): *Synechococcus* and *P.
    fluorescens* proteins containing one or two transmembrane domains
    according to TMHHM analysis of the whole genomes will be extracted
    and further analyzed using Psort and Pfam programs in order to
    identify and characterize a set of excreted proteins in these
    organisms.

-   Flagella and pili structural genes and related regulatory genes

-   Adhesion proteins

The results of identification and characterization of proteins relevant
to a process of biofilm formation will be passed to transcriptome,
proteome, and protein-protein interaction analysis groups. Information
regarding architecture of gene networks for these gene products and
composition of protein complexes under variety of experimental
conditions will be incorporated into the computational model of a
transition state between planktonic and attached fases of *Synechoccus
WH8102* and *P. fluorescens* existence.

**\
Deliverables**

 

***Year 1:***

-   Compilation and debugging of selected flux models based on the
    metabolic reconstructions

-   Fitting of the flux model of *Synechococcus* to different metabolic
    states. This involves measuring metabolic fluxes toward main storage
    compounds and main metabolic fluxes crossing plasma membranes.

-   Fitting of the flux model of *Pseudomonas fluorescens* to
    experimental data obtained at different modes of btach growth and
    physiology state

***Year 2:***

-   Sensitivity analysis of fitted models

-   Conversion of flux models into steady-state kinetic models;
    qualitative and quantitative analysis of kinetic models to revel
    defects in functioning uncontrolled central metabolism

-   Simulation of synchronous modes in homogeneous cell populations

***Year 3:***

-   Generation and debugging of stoichiometric matrices generated from
    the pathway collection, cross-linking with small molecules
    nomenclature

-   Simulation of dynamic behavior in homogeneous film of *Pseudomonas
    fluorescens*

-   Conversion of steady-state models of into dynamic models

***Year 4:***

-   Optimizing the dynamic models to predict control mechanisms
    governing intermediary metabolism of *Synechococcus* and
    *Pseudomonas fluorescens*.

-   Fitting the dynamic model of *Synechococcus* to experimental data
    obtained in experiments with continuous synchronous culture of this
    organism.

-   Conversion of the fitted dynamic model of *Synechococcus* sp. into a
    model describing of spatially homogeneous population of this
    organism

***Year 5:***

-   Simulation of dynamic behavior of a mixed homogeneous population of
    *Synechococcus* *and Pseudomonas fluorescens*

-   Simulation of spatially inhomogeneous cell populations
    *Synechococcus* and *Pseudomonas*

-   Simulation of dynamic behavior in homogeneous film of *Pseudomonas
    fluorescens*

#### References

Aponin 1971

\[Atkinson, 1969\] Atkinson, D. E. Regulation of enzyme function. Ann.
Rev. Microbiol. 66 (1969) 47-68

Avseenko, 1980

Avseenko 1987

Bioteux 1980

\[Boiteux, 1984\] Boiteux A., Hess B., Sel'kov E. E. -- Experimental and
theoretical studies of the regulatory hierarchy in glycolysis. In:
Self-Organization. Autowaves and Structures Far From Equilibrium. (V. I.
Krinsky, Ed.) Springer Series in Synergetics 28. Springer, Berlin, 1984,
pp. 240-251

\[Bunning, 1972\]

Cornish-Bowden 1995

\[Cowert, 2001\] Covert MW, Schilling CH, Palsson B. Regulation of gene
expression in flux balance models of metabolism. J Theor Biol. 2001 Nov
7;213(1):73-88.

Dunlap 1980

Ebeling 1983

Ebeling 1985

\[Edmunds, 1984\]

\[Edmunds, 1988\]

Edwards 2000

Enzyme Nomenclature 1992

\[Goldbetter, 1966\] or 1996?

\[Gurel, 1983\]

Heirich 1996

Hunding 1974

\[Ivanitsky, 1978\] Ivanitsky, G. R., Krinsky, V. I., and Sel\'kov, E.
E., Mathematical Biophysics of the Cell, Nauka Publishers, Moscow, 1978,
pp. 310 (in Russian)

Koshland 1966

Malkova 1980

\[Mano 1968\]

Mano, 1969

Mano, 1970

Mano 1971-76

Miakishev 1983

\[Michaelis 1913\]

Mood 1963

Morales 1967

Nazarenko 1983

Othmer 1976

Popova 1975

1976

1978

1979

\[Rapkine, 1931\]

\[Reich, 1975\]

\[Reich, 1981\] Reich, J. G., and Sel\'kov, E. E. Energy Metabolism of
the Cell. A Theoretical Treatise, Academic Press, London, 1981

Rensing 1980

Schneegurt 1994

Schulmeister 1978

Selkov 1967

\[Selkov, 1970\]

Selkov 1972

\[Selkov 1975\]

\[Selkov 1979\]

\[Selkov 1979a

Selkov 1979c

\[Selkov, 1980\] Kinetics of polyenzyme systems, DSc.Thesis, Russian
Academy of Sciences, Pushchino, 1980

Selkov 1980a

\[Selkov, 1980c\]

\[Selkov, 1980c4\]

Selkov 1989

\[Thorburn ,1915\] W. M., \"Occam\'s razor,\" *Mind*, 24, (1915), pp.
287-288

\[Tikhonov, 1952\] Tikhonov A.N. Systems of ordinary differential
equations containing small parameters for the derivatives. Mat. Sbornik.
31 (1952) 575-586

\[Vasilieva, 1967 \] A B Vasilieva, A.B. and Volosov, V.M. The work of
Tikhonov and his pupils on ordinary differential equations containing a
small parameter, *Russian Math. Surveys* **22** (2) (1967), 124-142