Volume 8, number 2


August 31, 1975

MATHEMATICAL MODELING OF DIFFERENTIATION IN DICTYOSTELIUM DISCOIDEUM David J. M. PARK and Barbara E. W R I G H T Department of Developmental Biology, Boston Biomedical Research Institute, Boston, Massachusetts 02114

(Received August 30, 1974)

Summary Methods for the dynamic analysis of biochemical differentiation are presented. These are demonstrated in the analysis of biochemical differentiation of the carbohydrate system in D. discoideum. Procedures for simplification which are presented are projection and contraction of the system trajectory in state space and the generation of reduced equivalent dynamic metabolic networks. The importance of the hierarchical structure of differentiating systems is discussed and the concept of a dynamic embedding diagram is introduced. It is shown that complex systems must be analyzed on an epoch by epoch basis, each epoch being a period of time characterized by a constant dynamic embedding diagram, and that widely different time scales and state space scales may be necessary in different epochs. In particular there is no a priori lower limit to the time scale which may be necessary during the analysis. Some problems in mathematically defining differentiation are discussed.

Introduction This paper explains some of the mathematical modeling techniques which the authors have developed for the study of differentiating systems. As these techniques are significantly new and different from many of the standard modeling techniques, their presentation appears to be justified. Although the general procedures used will be discussed, we will not attempt to present the details in this paper. Much of the formalism and computational procedures for modeling

biochemical networks are discussed by one of the authors elsewhere 1. In our laboratory these techniques are being used in the study of biochemical differentiation of the carbohydrate system in Dictyostelium discoideum. A complete description of the specific model used for this system is given by WRIGHT and GUSTAFSON2. For the purposes of the present paper, a simplified version of this system is used to illustrate some of the techniques. A paper published elsewhere 3 gives more complete calculations on the transformation of the D. discoideum carbohydrate system as it is affected by the various parameters of the system. Biochemical systems are of such complexity that one is eventually forced to resort to the language of mathematics to describe them. A set of variables and the mathematical equations which relate them constitute a mathematical model for a system. The variables in the model must be placed in a correspondence with the real physical and chemical entities in the system. The better the mathematical model reproduces the real behavior of the system, under a variety of conditions, the better the model is considered to be. The mathematics which deals with the dynamic interrelationships of a reasonably large number of variables is called systems theory. The language of systems theory can be used or we can invent our own equivalent chemical language. We regard the first option as preferable. In this regard much of the standard terminology is explained by ROSEN4. There is one rule which must be followed. Any concepts entertained must have their correlation in the mathematical model. A prescription must

Dr. W. Junk b.v. Publishers - The Hague, The Netherlands


be available which describes the method of calculating a quantity from the model and, conversely, the model should suggest how the quantity could be measured in the actual system. Concepts and quantities which do not have a representation in the model should be abandoned, or else the model should be extended and prescriptions given for obtaining them from a model. Quantities obtained from a model are, of course, only an approximation to the real quantities. A more complete model or more accurate calculation may alter their values, but, hopefully, this is a convergent process. Many of the techniques are presented in terms of a geometrical representation. Hence, the systems may be thought of as points which move along paths or trajectories in a high dimensional state space. These trajectories describe the changes in chemical composition and configuration of the system as a function of time during differentiation. With this geometrical picture in mind, questions such as the following may be asked: When a parameter of the model is changed, how does it affect the trajectory in state space? Does it affect the direction of the trajectory? Does it affect the speed with which the system moves along the trajectory? How much will the system deviate from its normal end point? Since any interesting biological system is extremely complex, techniques of simplification must be invoked in order to make analytical progress. There are valid techniques for simplification which allow a system to be transformed into a manageable state from a mathematical standpoint. The methods to be presented, however, are different from those currently in use, such as altering rate constants to increase the stability of the solution procedure or arbitrarily selecting a subsystem for study by its chemical as opposed to its dynamical character. Techniques to be discussed include those of projection into substate spaces accompanied by contraction of similar substate spaces to produce bundles of trajectories, and metabolic network reduction based upon the detection of steady state subnetworks (which are suppressed), conserved moieties (which are parameterized), and dynamic moiety pools (which become new variables). This paper will concentrate upon the dynamical character of differentiating systems. Such systems have a hierarchical structure and this structure 98

must be taken into account in order to work with efficiency and to discover the dynamically important variables during each epoch of a process. This will lead to the concept of the active hierarchical level of a network and allow the identification of variables associated with this level as the critical variables of the network during the epoch. This concept has already been presented in biochemical terms by one of the authors 5, but is presented here in the context of the mathematical theory of hierarchical systems. Usually, in solving differential equations, it is natural to think in terms of uniform time and space. A time span of interest is chosen and divided uniformly; a similar process is applied to the range of chemical and other variables. It should not be surprising that this procedure is inappropriate for biological and biochemical systems. Advanced simulation procedures will allow extremely nonuniform time and space scales. The fineness of these scales is a function of the hierarchical level and is dictated by the dynamics rather than being freely open to choice. Lastly, the difficulty of mathematically distinguishing °differentiation' from other types of biochemical and biological transformations will be discussed. Sometimes it appears that the very process of modeling a 'differentiating' system robs it of its essential characteristics. When these processes are modeled they begin to look like transients and it may be difficult to believe that 'differentiation' is a mere transient. But this is because transients are usually considered in systems with a narrow hierarchical range. In highly nonlinear systems with a wide range of hierarchical levels transients begin to look much more like the complex processes commonly associated with 'differentiation'. As a complex process proceeds from one epoch to another the set of variables which comprise the active hierarchical level will change. Thus, different variables become dynamically active at different times, and it should not be surprising that the process cannot be described by selecting a single class of variables, say genetic states. The hierarchical structure of metabolic networks tends to shield information stored at lower levels from the active hierarchical level. Thus, as far as determining the behavior of a system during any given epoch is concerned, the amount of information stored in the dynamic variables may

be comparable to the information transmitted, for example, from the genetic variables.

Geometrical Representation Geometrical representations are very convenient for the discussion of the behavior of complex dynamical systems. To proceed with such representation a state space must first be constructed. A state space will have one axis for each variable in the system. The state space can be visualized as completely orthogonal and Euclidean. The state of the system can then be represented by marking off on each axis the value of the variable which corresponds to that axis and then plotting the resulting point. Hence at any given time the system is represented by a single point in state space. However this point will generally move as a function of time and hence a process can be represented by a trajectory in state space on which we can place tick marks to represent the time intervals along the trajectory. To put some substance on this picture, consider the D. discoideum system and the construction of a single fruiting body. Beginning with an aggregation field of, say, 10,000 cells, the xyz spatial coordinates of these cells can be listed as a function of time. With just these three variables for each cell there are already 30,000 variables and hence the state space is already of quite high dimension. Next the chemical composition of each cell must be described but before that is done it is necessary to recognize that each cell may be divided into physical compartments and phases. As an example, let us say that each cell is divided into 100 compartments and we use 100 chemical variables to describe each compartment. The chemical variables would be the pool sizes (number of molecules) of each metabolite in the compartment or phase. This now gives approximately 108 variables to describe the system. There may be additional variables such as those associated with spatial diffusion fields of external metabolites. The large number of variables associated with such a description presents severe problems. Without simplification, the prospects of getting anywhere with modeling are rather dismal (or with any other kind of investigation; if it is difficult to simulate 108 variables, imagine measuring all of them !). As it is difficult to visualize a high dimensional state space, visualize it in two or three

dimensions. The mathematics is devised to handle any number of dimensions and is a purely mechanical process. The chemical state space for each cell is the union of the individual chemical state spaces of each compartment in the cell. The chemical state space for the entire system is the union of the chemical state spaces of each cell. The state space for the entire system is the union of the chemical state space and the xyz spatial coordinates of each cell. A complete mathematical description of a differentiating system would consist of a set of simultaneous nonlinear differential equations, one equation for each variable, and these equations would contain only the variables themselves and constants. The constants which appear in these equations can be called the parameters of the system. The concept of state space can now be extended to include an axis, not only for each variable, but also for each parameter. Questions may then be asked concerning the effect of the parameters on the outcome. If one parameter were different this would correspond to starting the trajectory from a slightly different point in state space. This would then give a different trajectory and a different outcome. In incomplete models, such as those to be presented, the parameters are not all constants but some are given by time functions. This will then give parameter trajectories in the parameter portion of state space. Many of our analytic studies have been concerned with determining the system trajectory under the action of various parameter trajectories.

Projection and Contraction A specific instance of differentiation in a specific system, say the construction of a single fruiting body from the aggregation field of amoeba in D. discoideum, is represented by a single trajectory in a high dimensional state space. This is illustrated schematically in level (A) of Figure 1. The trajectory starts at some point A (vegetative amoeba), snakes and twists its way through 108 dimensions of state space, and terminates at a stable steady state position B (fruiting body). A time scale is placed upon the trajectory by the use of tick marks~ The point A may have initially been a stable steady state which is then made unstable by an exterior event (deprivation of food supply in D. discoideum). The first method for simplification of such a 99



L:> ....... E......

CELL I0,000









Fig. 1. Projection of complete system trajectory. Only three of many dimensional axes are shown in each plot. The axes represent the variables of the system but are left unlabeled because there are many others not shown. (A) Complete system trajectory. (B) Projection into the chemical state spaces of each cell. Only three of the 10,000 cells are shown. (C) Projection from one cell into the chemical state spaces of each compartment in the cell. Only three of the 100 compartments are shown. (D) Projection from one compartment into the state spaces of subnetworks. complex system is simply to select out certain variables for study. Geometrically this corresponds to the projection of the system trajectory into a smaller dimensional state space. This is similar to the projection of a three dimensional curve into a two dimensional plane. When this is done some information is lost (movement in the third dimension), but the information could be preserved by doing a series of projections. Returning to Figure 1 we can project the total system trajectory into the chemical state spaces of each of the 10,000 cells. This is shown schematically in level (B). Each one of these state spaces m a y have approximately 104 dimensions. Next the trajectory of each cell in its own state space can be projected into the chemical substate spaces which describe the compartments of the cells, This is shown in level (C) of Figure 1. Each 100

o f these substate spaces m a y have, say, 100 dimensions. Lastly, the chemical variables of a specific c o m p a r t m e n t m a y be further projected into the state spaces of specific subnetworks (D), for example, the network interconnecting the major carbohydrates in D. discoideum. Projection is a necessary process in arriving at a manageable subsystem for study. However, it is also a dangerous process if it is done arbitrarily and without regard to dynamic considerations. I m p o r t a n t variables may be thrown out and never again considered, or variables which are tightly linked dynamically m a y be separated. Projection can, however, be a very useful process if it is combined with another geometrical operation called contraction. This is illustrated in Figure 2, and starts with a large number of subsystems which have been projected out of the total system trajectory. These subsystems might be a large group of cells from level (B) of Figure 1, or they m a y be a number of similar compartments, say mitochondria, from level (C) of Figure 1. They m a y be a number of similar subnetworks, say carbohydrate metabolism, from level (D) which have been projected out of each o f their cells. In each case the subsystems contain a similar n u m b e r of dimensions and similar variables. These similar state spaces can now be superimposed upon each other making sure that




....... r




IP /



Fig. 2. Contraction. (A) A number of similar subsystemsfrom level (A), (B), or (C) of Figure 1. (B) Superposition of similar state spaces when trajectories are identical (C) Superposition of similar state spaces when trajectories are similar but not identical.

similar axes are aligned (B). This process we call contraction. When this is done it may turn out that all of the subtrajectories lie on top of each other as shown in Figure 2(B). When we say that the trajectories lie on top of each other, we also mean that the equal time tick marks also lie on top of each other. If a picture such as shown in (B) occurs we have then succeeded in collapsing the system state space into many fewer dimensions, cashing in m variables for 1 variable, and doing this for each of the chemical components of the subsystem. Of course, superposition of all of the individual trajectories may not result in a picture that looks like Figure 2(B). The trajectories may not lie on top of each other, but may be distinct. They may form a bundle of trajectories that looks like the picture in Figure 2(C). Each subsystem (cell, organelle, subnetwork) follows its own trajectory. The equal time tick marks can be connected into a set of surfaces which cut across the bundle of trajectories. Call these the equal time surfaces of the bundle. If the bundle is too fat, or if the equal time surfaces are too skewed (i.e., do not cut across the bundle at right angles) then the bundle must be divided into more than one group. Divide a fat skewed bundle into a number of smaller bundles which are not fat or skewed. Figure 3(A) shows a well behaved bundle of trajectories. The figure contains a discrimination box indicating how finely the chemical composition is to be specified. The box fits over the bundle of trajectories so they can be regarded as all the same within the accuracy requirements. Figure 3(A) also shows the equal time surfaces which cut across the bundle of trajectories. They are orthogonal to the axis of the bundle. This is the way it should be. Figure 3(B) shows a bundle which is properly thin but which is skewed because the equal time surface slices the bundle at an acute angle. This means that although all subsystems in the bundle are moving along a similar path, they are moving along it at different rates and times. This does not form an acceptable bundle (some subsystems may 'miss the bus') and must be further subdivided. Figure 3(C) makes this concept more precise. On the left a bundle fits within the discrimination box. The equal time surface cuts the bundle at a slight angle yet the extreme points lie within the discrimination box. This would be acceptable. On the right, the bundle again lies within the discrimination box.











> Fig. 3. Skewness of trajectory bundles. (A) A thin nonskewed bundle. (B) A skewed bundle. (C) Relation of discrimination boxes to skewness.

However, the equal time surface cuts the bundle at such a skewed angle that the extreme points of cutting lie well outside of the discrimination box. This would not be acceptable. It may appear that the process of dividing up bundles is one that can be carried out intuitively and by inspection. This, however, is wrong as will be discussed below. The size and shapes of the discrimination boxes are not arbitrarily open to our discretion and can be a strong function of their location in state space. A complete description of a differentiating system contains both the three dimensional movements of the cells and the movement in chemical state space. This is illustrated in Figure 4 for D. discoideum. If the single system trajectory is projected into the xyz spatial coordinates of each cell and then these physical state spaces are contracted the picture of Figure 4(B) is obtained, something like a smeared motion picture of the 101















Fig. 4. A completedescription requires both morphological and chemical variables. (A) Singlecompletesystem trajectory. (B) Picture obtained by projectingthe completesystem trajectory onto the spatial coordinates of each cell and then contracting these spatial coordinates.This is the usual picture of the morphologicaldevelopmentof a single fruiting body from an aggregationfield of amoeba in D. discoideum. (C) Picture obtained by projectingthe completesystem trajectory onto the chemical coordinates of each cell and then contracting these chemical coordinates. development of a single fruiting body. However, if the single trajectory is projected into the chemical state spaces of each cell and then contracted a picture such as Figure 4(C) is obtained. Here all amoeba start out pretty much the same at point A but form two distinct cell types, spore cells and stalk cells. The stalk cells themselves may be chemically graded and are therefore shown as spread out. Motion pictures have been made of Figure 4(B)* and are well known to workers in the field. No motion picture has ever been made of Figure 4(C), although much more information is contained here, and this is just as much a part of the differentiation process. One way to 'explore' state space is to calculate a large number of trajectories for slightly different starting positions. Biology is reduced to ray Princeton University,Princeton, N.J. and COHEN,M. H., Universityof Chicago, Chicago, Ill.



tracing. It is much easier to do this on the computer than in the laboratory, especially if we just want to check out an idea. State space does need 'exploration'. There is no way to solve a large number of highly nonlinear simultaneous differential equations by intuition or inspection. Simply having a model, even if it is correct, does not mean that its behaviour is understood. Thinking again of the trajectory of the system in the complete, unreduced, state space it is possible to visualize the process as similar to a shooting match. The final position can be thought of as a 'target' which the system must 'hit' for successful differentiation. In fact, thinking in terms of the contracted chemical state spaces of the cells (Figure 4(C)) there is clearly not one target but many - and the system succeeds in hitting them all ! Nature learned how to MIRV its cells long before we learned how to MIRV our rockets. Researchers often talk about the 'control' of differentiation. It is meaningful to talk of control only it if is clearly stated what is being controlled and with what accuracy. In D. discoideum this might be expressed in terms of how closely the final composition and structure of the spore and stalk cells is controlled. Perhaps the spores must be more closely regulated than the stalk cells. Certain chemical variables may be much more closely regulated than others for the system to be viable. Unfortunately, this is an area in which there is very little quantitative information. Suppose, however, that this information was available. There might be 'surfaces' in state space which represent constant probabilities of survival for the system. This would then give some indication of the size of the discrimination box needed at the end point of the system trajectory. Knowing this, the discrimination box could be slid back along the system trajectory, but as this occurred its size and shape would have to change in such a way that any trajectory which passed through the box at an earlier time is assured of ending in the final box. This is illustrated in Figure 5. The accuracy with which the system must be characterized is not freely open to choice, but is dictated by the dynamics of the process. Of course, it is not necessary to explore all of state space, but only those regions near to the actual trajectory. It would be desirable to have the trajectory itself and some complex quantities



1532" A

Fig. 5. Discrimination boxes are a function of position on the ;ystem trajectory.

which go along with the trajectory. One of these quantities would be the size and shape of the discrimination boxes (to hit the target) at each point on the trajectory. Another important quantity associated with the trajectory at each point is the gradient of the net flux of the metabolic pools (or the gradient of the rate of change of other variables). This gives information about stretching along the fiducial trajectory, the skewness of nearby trajectories and even twisting of the nearby trajectories about the fiducial trajectory. Since advance simulators will often need many of the components of the gradient, this is information which should be easily available. Metabolic Network Reduction

By collapsing state space and forming bundles of trajectories it is possible to reduce the complexity of a system and the number of independent equations which must be solved. Nevertheless, the number of metabolic variables within a single cell or compartment may still be so large that severe problems remain. Fortunately, there are further possibilities for simplification. There is a general procedure published elsewhere 1 giving the conditions under which a metabolic network can be reduced to a simpler network. The procedure is based upon the steady state approximation. It utilizes an error criterion and generates a set of reduced equivalent

networks. The net result of the procedure is to suppress some pools from the network, and to form new pools from a linear integral combination of the old pools. There is a net overall reduction in the number of dynamic variables. This process will be briefly described and then illustrated on several networks. A steady state pool is essentially one in which the net flux (inward minus outward) is much smaller than the outward flux. The steady state condition is a dynamic relationship between a metabolic pool (or subnetwork) and its embedding network. This relationship allows the steady state pools to be solved algebraically by setting their net flux equal to zero and solving the resultant set of nonlinear algebraic equations. These algebraic equations replace the original differential equations representing the pools. These differential equations happen to be those that cause great difficulties in performing efficient simulations. They represent what are known as 'stiff' differential equations 6, and are a serious pitfall for anyone contemplating metabolic simulation. Whenever possible it is well to be rid of them. Solving a set of nonlinear simultaneous equations is not easy; for some guidance see Traub 7. The process is essentially one of guessing and improving by an iterative procedure. If the initial guess is too far off there can be real computational problems. However, in a real simulation process the picture is much brighter. When the pools are not in steady state, it is not necessary to solve the algebraic equations; the differential equations are solved instead. When the pools do approach the steady state condition and it is necessary to switch over to the solution of the nonlinear algebraic equations, a good approximation to the solution is already at hand. Steady state subnetworks can simply be removed from the dynamic network. When this is done some reaction lines connected to the remaining pools are left dangling. How they are reconnected depends first on how many independent reactions remain in the reduced network, and second on which subset of the original reactions we choose to build the reduced network. This means that it is possible to arrive at a number of equivalent reduced networks. A set of connected steady state pools may contain within themselves nonsteady state dynamic variables. These variables are linear positive integral combinations of the steady state 103

pools and are defined as moiety pools. A moiety pool may be absolutely conserved within a steady state subnetwork in which case the pool may be parameterized, that is, it may be treated as a constant parameter of the system. More generally, the moiety pool is a dynamic pool which must be added back to the set of differential equations. Contrast this mathematical procedure with one much more familiar to biochemists, the derivation of enzyme rate equations. There, one starts with a set of linear equations to solve, which is much easier. Perhaps because of this the use of the steady state approximation is quite c o m m o n at the enzyme level, and yet there appears to be a psychological barrier to extending it above that level, despite its complete justification on dynamical grounds. In the case of the derivation of enzyme rate laws, the enzyme moiety plays the role of a conserved moiety, although in m a n y cases it would more accurately be represented as a dynamic moiety pool. As a first illustration, consider the network of Figure 6. This shows two metabolic pools which are in rapid 'equilibrium'. Imagine that the reactions involving the two pools are well approximated by mass action rate constants as shown. Then the two pools will reach 'equilibrium' in a matter of milliseconds. If the input flux, go, which represents the effect of the embedding network, varies slowly enough then the response of network (A) will be well approximated by the network (B). However, if the input flux contains a sharp discontinuity then network (A) must be used. The network reduction rests purely on dynamic grounds. In the network reduction of Figure 6 the two steady l k : 1000 mTn-1 k 3 - 1 mM-1

k 2 = 1000 min -I (A) go





= 0.5



Fig. 6. A steady state subnetwork with a dynamic moiety pool. The network is driven by the input flux go. The metabolite concentrations are given by xl and x2. The mass action rate constants are given by k~, k2 and k3. (A) The dynamic network representation. (B) The reduced dynamic equivalent network when the two pools are in steady state or 'equilibrium'. The two original pools disappear and are replaced by a single moiety pool with concentration xa + x2. 104

' L ~



'q R4 R4


3 OR



R R5










Fig. 7. (A):A simplifiedrepresentation of carbohydrate metabolism in D. discoideum between culmination and sorocarp. The metabolic pools are labeled P1 through P6 and the reactions are labeled R1 through Rs. The biochemical identification of these pools and reactions is given in Table 1. During this period the pools in the dotted box forin a steady state subnetwork with no conserved or dynamic moiety pools. (B) to (G) : Reduced equivalent dynamic networks of (A). A number of equivalent networks arise because there is choice in which reactions of (A) are to be used in the reduced network. If a reaction involvesa stoichiometric constant other than flnity then the constant is shown near the end of the reaction line. In (E) the P3 pool is weighted by a factor of two to avoid fractional stoichiometric constants. state pools are suppressed from the network but a new moiety pool is introduced as a dynamic variable. Figure 7(A) presents a simplified network for a major portion of carbohydrate metabolism in D. discoideum cells between the culmination and sorocarp stages of differentiation. The network is an approximation, but not a bad approximation, and allows some of the methods of analysis to be illustrated without wandering in state spaces of too great a dimensionality. In the network the pools are labeled P1 to P6 and the reactions are labeled R 1 to R s. The biochemical identification of these metabolic pools and reactions is given in Table 1 along with some pertinent information

Table 1 Metabolic pools and reactions of the network of Figure 7(a) Metabolic Pools: Px is soluble glycogen measured in terms of the number of glucose units. P2 is cellulose wall material also measured in terms of the number of glucose units. Pa is trehalose P4 is glucose-l-P (G1P) Ps is glucose-6-P (G6P) P6 is uridine diphosphate glucose (UDPG) Metabolic Reactions: R1 is glycogen phosphorylase (EC Mechanism is rapid equilibrium random bi bi with B=Q K(inorganic phosphate) = 3 mra K(sol. glycogen) = 3.2 mM K(G1P) = !.2 mM V1 given by time function Keq = 0.33 (Wright & Gustafson, 1972) Ki(G1P ) = 0.1 mM R2 is cellulose synthetase mechanism is linear k = 0.52min l R 3 is trehalose synthetase (EC mechanism is ordered bi bi with P, Q ~ 0 K(G6P) = 10 mM K ( U D P G ) = 2.2 mM Ki(G6P) = > 50 mM V~ given by time function (Wright & Gustafson, 1972) R 4 is UDP-glucose pyrophosphorylase (EC mechanism is ordered bi bi with P ~ 0 K(uridine triphosphate) = 0.1 mM K(G1P) = 0.5 mM K~(uridine triphosphate) = 0.1 mM Ki(UDGP) = 0.05 mr~ V~ = 45.9 mM/min R 5 is phosphoglucomutase (EC in rapid equilibrium, Keq = 6.0

concerning the reaction mechanisms employed in the model. Under the conditions that actually occur in this network the set of pools (P4, P 5 , P6) enclosed in the dotted box form a steady state subnetwork. Under the physiological conditions that occur, the use of the steady state approximation does not alter the calculated values of the pool sizes and reaction rates by more than 1 or 2 ~ . Similarly, these steady state pools never contain more than 2 to 3 ~ of the total mass of glucose moieties in the network. This can either be ignored or the material can be absorbed back into the embedding network before carrying out the network reduction. In this steady state subnetwork there are no conserved or dynamic moiety pools (as there were in Figure 6).

The exact numerical procedure for generating the reduced equivalent networks for (A) or any other network is given elsewhere 1. Here only the results are presented. It may be noted, however, that the procedure is neither trivial nor especially intuitive. For example, jump immediately to Figure 7(F) and (G) and see if those networks appear to be reduced networks of (A) or if it is obvious that (F) and (G) are completely equivalent to each other. Removal of the steady state subnetwork from (A) leaves some dangling reaction lines which must be reconnected. In this particular case there are two independent reactions in the reduced dynamic network, so when the remaining pools are reconnected there will be two reactions. There are 10 ways of picking two reactions out of a set of five, suggesting that there are 10 possible reduced equivalent networks to (A). However, under the steady state condition, R 3 and R 5 are linearly dependent reactions and cannot be used as a set, so nine different reduced equivalent networks remain. Reactions R 3 and R 5 always enter into the reduced networks in the same way so there are actually only six different forms of the reduced networks and these are shown in Figure 7(B) to (G). It should be stressed that the reaction rates to be associated with the reduced equivalent networks are the same as those which occur in (A). For example, the rate of R4 in (C), (F) and (G) is exactly the same as the rate of R4 in (A). In these networks the stoichiometric coefficients are shown as integers near the end of the reaction lines unless they are unity, in which case they are suppressed. In network (E) the P3 pool is weighted by a factor of two to avoid fractional stoichiometric coefficients. The second part of the procedure for reducing networks consists of the actual solution of the steady state equations. This solution plays the same role as the 'rate equation' does in enzyme kinetics. The actual derivation of the solution for this network is given in the Appendix and the result of that calculation is given in Table 2. Here the rate of reaction R1, r 1, is regarded as an independent variable and all other reaction rates are displayed as functions of r 1. Note that the steady state subnetwork acts both as a flux divider (r 1 = r 2 + 2r3) and also as a flux limiter t o P2 (cellulose) because at high input flux levels the flow to P2 tends to saturate and all excess flow is diverted to P3 (trehalose). 105

Table 2 Calculated Response of Flux Divider Network* rl 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000

r2 0.0791 0.1234 0.1509 0.1697 0.1836 0.1943 0.2028 0.2098 0.2157 0.2207

r3 = r5 0.0104 0.0383 0.0745 0.1151 0.1582 0.2028 0.2486 0.2951 0.3421 0.3897

r4 0.0896 0.1617 0.2255 0.2849 0.3418 0.3972 0.4514 0.5049 0.5579 0.6103

*k0 = 0.13 mM; k~ = 1.05 mM/min; k 2 = 2.85 mM- 1 m i n - i ; k3 = 0.52 m i n - i ; r~ through re in mM/min.

To see that these networks are all actually equivalent one can calculate the net flux on each of the pools in all of the networks and see that the same answer is always obtained. For example, the net flux on P~ is given by the following expressions in the various networks: (B) - - r 2 -- 2 r 3 , ( C ) - r 3 - re, (D) - r , , (E) - r x , (F) r2 - 2r4, (G) - r 1. By taking a particular set of values from Table 2, for example r 1 = 0.5, r 2 -- 0.1836, r 3 = 0.1582 and r 4 = 0.3418 and substituting into the above expressions it can be seen that the net flux on Pa is - 0 . 5 for all networks. The same equivalency is obtained for all pools and all sets of values in Table 2. Therefore these networks are all equivalent because equivalence is defined as the same net flux on the dynamic pools. These networks are not appropriate for answering all biochemical questions; only (B) represents a correct precursor-product relationship for tracer experiments. But the other networks can give insight into the network behavior not given by (B). Network (G) at first sight looks preposterous. Pa flows to P2 ? and via Re !? Nevertheless (G) is not only perfectly correct but also useful. Suppose the Rj reaction is held fixed and the Re reaction is activated. What will happen? Network (G) predicts that there will be an increase in the net flux o n P 2 and a decrease of the net flux o n P 3 . P 2 will build up at the expense of P3 (not at the expense of P~). Figure 8 shows that this is actually the case. Here the results of the reduced metabolic network are combined with representation of a system response as trajectories in state space. The system 106

trajectory has been projected into the P2 - P3 (cellulose-trehalose) substate space. These trajectories were actually calculated without using the steady state approximation. Also the R 1 reaction undergoes a time dependent deactivation during the calculation but this is exactly the same for all trajectories. The different trajectories represent different activations of R 4 and show that indeed material is shifted from trehalose to cellulose as R 4 is activated. The state space plot of Figure 8 is, incidentally, nearly a complete representation of the system response. Because of the conservation relation on the glucose moieties it is not necessary to plot the concentration of P1 (glycogen) as that would add no new information. In Figure 8 the trehalose concentration is weighted by a factor of two in order to be expressed in terms of glucose moieties. In this case the 45 ° dashed lines represent surfaces of constant concentration of P2 + 2P3- Because of the conservation relation they also represent surfaces of constant P~ concentration. However, they were actually drawn in to give the constant

300 MI N

R4/40OR 40R3

25oA~', \



\,>% ~.


75/. -, ,..,,.... ,J. "\ \\

~ ~_~

I ' \ ".,. \ \ \ " " , ~51 \ \ \ \ \ ~

z z









,5 °4. " . \ . \ . Y 2 \ % ,o

\\\ / /\ \ \\\,\ \ \. \

4, I



o 5

\\\'\R 4'R.~














Fig. 8. Carbohydrate trajectories in D. discoideum. This is a projection of the trajectory of the network of Figure 7(A) into the cellulose-trehalose subspace. The concentrations of cellulose and trehalose are both given in terms of glucose units and on a packed cell volume (p.c.v.) basis. The dashed lines represent equal time surfaces for the trajectories and are calibrated in minutes after culmination. The fiducial trajectory is the one most nearly agreeing with the experimentally obtained concentrations. The other trajectories represent ones obtained by activations or deactivations of reactions R 3 a n d R4. A similar result can be obtained by either activating R a or deactivating Re by the same factor,

time surfaces of the trajectories. This shows that indeed the extra accumulation of cellulose is at the expense of trehalose and not glycogen just as network (G) predicted. The amount of glycogen used in a given amount of time is the same for any of the trajectories. There is a slight departure from this on the R4/40 trajectory at early times because there the steady state approximation is starting to break down and a significant amount of material is building up in the steady state pools. The other networks are useful for answering other questions. For example, Table 2 shows that the rate of reaction R2, r2, increases more slowly than the rate of R 1. What happens if the rate of glycogen breakdown is increased and r 2 is held constant? Network (E) gives the answer to this question. All the extra flux goes to trehalose. Or what happens if R 1 is held constant and R a is deactivated? Network (D) gives the answer to this question. P2 is increased at the expense of P3" This is the same effect as occurred in (G) when R 4 was activated. These two effects are not only qualitatively similar but the Appendix shows that they are quantitatively equal. It is only the ratio of the R 4 and R a activities which is important. Hence Figure 8 is also a good representation of this parameter variation. Another comment should be made about the steady state approximation and functions like those shown in Table 2. From Table 1 it can be seen that there were a fair number of kinetic constants that went into the generation of the functions of Table 2. Because of the equivalency relations only r 2 is needed as a function of r 1. This r 2 function could be given to a mathematician with the request: 'Approximate this curve over the range shown to an accuracy specified by using an easily evaluated equation with as few constants as possible.' We might be pleasantly surprised to learn that the result contains many fewer constants than there were kinetic parameters that went into the original generation of the function. This investment could then produce a model which is simpler to handle. As important as this is for simplifying the calculations, it is even more important to the experimentalist. Considering the cost of biochemical experimentation, it is highly desirable to obtain some idea of the number of important dynamic parameters needed to characterize a network.

Dynamical Characterization of Differentiating Systems The introduction of the steady state approximation, and consideration of the problems in solving 'stiff' differential equations gives a new way to characterize the various metabolic pools in a metabolic network. The characterization is on dynamical grounds. First there are the steady state pools. These pools drop out of the dynamic representation of the network because they are locked to the dynamic pools which embed them. The steady state pools can be calculated algebraically from the values of the dynamic pools. The steady state pools have the potential for rapid response but they are never challenged, and hence never have the opportunity to express this potential. Second there are the dynamic pools themselves. These pools are described by differential equations. They determine the timing of the metabolic transformations. The dynamic pools correspond rather closely to the concept of critical variable as one which is limiting the rate of differentiation during a particular epoch. Finally, there are pools which are so large and dynamically inert that they can be regarded as essentially constant over long periods. They are described fairly well by differential equations whose solutions are constants. So at any given time it is the dynamic pools which are interesting and important. These pools form a distinct subset of all metabolic pools in a given network and often they may be a very small subset. Hence again the modeling prospects are not as bad as they might seem without these considerations. The steady state pools are envisioned as associated with the low hierarchical levels in a metabolic system, the dynamic pools with the intermediate hierarchical levels, and the inert pools with the high hierarchical levels. This is in agreement with the concept that low levels of hierarchies are associated with short relaxation times and high levels with long relaxation times s. Hence a metabolic system can be thought of as having this threeway division. It should be noted that the steady state pools may actually encompass many hierarchical levels and the same is true of the inert pools. We call that portion of the system which consists of the dynamic pools, the active hierarchical level of the network. Figure 9 illustrates the division of a 107

Fig. 9. A dynamical embeddingstructure for a metabolic network. The diagram must be imagined as superimposedon a metabolicpathways chart of some complexityso that some pools lie in one region and some pools in another. The inert pools lie in the I region. The dynamicpools lie in the A regions. The steady state pools lie in the S regions. There may be multiple regions of each type and S regions may be directly embedded in I regions. The structure will evolveas a function of time; but will remain constant during any given epoch. hypothetical system into the three basic hierarchical levels. Imagine this diagram being superimposed on a metabolic pathways chart of some complexity so that some metabolic pools were situated in one region and some in other regions. The region containing the inert pools is denoted by 'I'. The active hierarchical level is denoted by 'A'. The steady state level is denoted by 'S'. This diagram shows that the active hierarchical levels are embedded in a constant environment and that the active hierarchical levels themselves embed the steady state levels. There may be several active networks which are disconnected from each other and each active network may embed none or many steady state subnetworks which are disconnected. We assert that this dynamic embedding structure of a network is just as important in characterizing differentiating systems as the more standard characterizations. The familiar metabolic charts give only a vague, if any, indication of this dynamic structure. To complicate this picture even further this dynamic structure must not be thought of as a constant feature. Indeed, in a complex process the dynamic structure will vary so that the active hierarchical level will range over many of the hierarchical levels of the network as a function of time. Associated with these variations in the dynamic structure will be variations in the time response of the network. When the active hierarchical level is at a high level the network 108

will be responding slowly. When the active level is low, the network, or parts of it, will be undergoing rapid transformations. This means that the time scale used in advanced simulations will be extremely nonuniform. At some epochs of the process it will be very fine, at other epochs very coarse. This time scale could vary by orders of magnitude. The output of advanced simulators will then look quite different from current simulators (including our own). Instead of displaying all variables on a uniform or nearly uniform time scale, only active variables will be displayed on a very nonuniform time scale. With these considerations in mind, the term 'epoch' can now be given a more precise definition as a period of time characterized by a constant dynamic embedding structure. It sometimes appears convenient to characterize systems by the time scale on which they operate. C. H. WAOOIN6a'ON9 divides biological processes into biochemical, epigenetic, and evolutionary and this division is followed by B. C. GOOOWIN1°. This division is dangerous if it implies that there is a certain lower level in time scale which can be ignored. There is not. For example, in dealing with a process which takes 24 hours, processes which occur on the time scale of milliseconds cannot necessarily be ignored. When the steady state approximation breaks down it must be abandoned with the concurrent acceptance of a finer time scale. There is no a priori lower limit. As a hypothetical example, consider again the network of Figure 6 and its steady state reduction. Imagine that one of the pools acts as an allosteric inhibitor and the other as an activator for an enzyme which catalizes an essentially irreversible reaction, so that the net reaction rate depends upon the ratio xl/x2 and only proceeds when xl/x2 is much greater than the equilibrium value. This reaction might be one which transforms an inactive enzyme into an active one. A sudden perturbation to the network of Figure 6(A) would effect this. Even if this all occurred very rapidly it would inevitably move the system to a new trajectory which in succeeding time would cause it to diverge more and more from its initial trajectory. If the steady state approximation of Figure 6(B) were always used, this process would never occur and an entirely different, and incorrect, answer would result. For this reason a division based upon epochs is preferable to one based on an arbitrary time scale.

Workers in gravitational physics use the phrase 'many fingered time' to indicate that it is often useful to push a solution forward at different rates of time in different regions of space 11. A similar procedure can be useful in the dynamic solution of metabolic networks. Looking again at Figure 9 it is seen that there can be several active networks which are disconnected. Remember also that Figure 9 represents a snapshot of the dynamic structure at a particular instant of time. This figure should really be a motion picture so that with passing time some of the 'A' pools may merge into the 'S' pools; some of the T pools may become 'A' pools; the two 'A' networks may later merge together, etc. Imagine studying a 24 hour process, and that the A1 and A 2 networks have just popped up as the active level. The A1 network may have a relaxation time of 1 minute and the A 2 network may have a relaxation time of 1 second. It would be very inefficient to push both the Aa and A 2 solutions forward at the same time rates. It would be much more convenient to use 'many fingered time'. This would be especially so if computers could be used which allowed multitasking computations.

Defining Differentiation Although we are investigating differentiation, the mathematical methods are applicable to a much wider class of phenomena than those usually characterized as differentiation. The methods apply to all types of biochemical, cellular, and biological transformations. From a mathematical and dynamical standpoint it is difficult to define differentiation in a manner which uniquely distinguishes it from other types of metabolic transformations. One possible criterion could be the number of variables in the system. For example, a single metabolic pool undergoing a transient between two steady states would not generally be considered a differentiating system. But how many variables must there be before it is a differentiating system? If some number such as 1000 variables were picked as the dividing line then no one is studying differentiation. The simulations with the largest number of metabolic variables are apparently those of GARFINKEL* * Personal communication.

which contain up to 300 metabolic pools. But he only claims to be studying intermediary metabolism ! Furthermore characterization by number of variables ignores the fact that there may be only a small number of dynamic variables at any given time. All metabolic transformations are nothing more than a transient from one steady state operation point to another. This should be stated because many studies on differentiation, especially those involving modeling, are often discounted as merely studies of transient behaviors. Mathematically, there is simply a set of simultaneous differential equations and their initial conditions. These equations are solved until a final stable state is reached. That this looks like something more than a transient can only be due to the extremely nonuniform time scale and state space scale which must be used. Another potential mathematical criterion which might be used to distinguish differentiating systems from transient systems is the following. If the differential equations all have constant parameters, then we might say only transient behavior is being studied, but if the equations have time dependent parameters, then the process involves differentiation. This accords with common concepts. For example, the enzyme maximum velocity parameters are dependent upon the amount of enzyme. If the latter changes because a gene is 'turned on' and new synthesis results, these parameters will be time dependent. Hence the 'essence' of differentiation may be defined in these time dependent parameters; they are what must be studied. This is a mathematical facade. The time dependent parameters appear only because some of the differential equations which describe the system have been suppressed. The time dependent parameters can be eliminated by writing the equations for the variables they represent. For example, it is possible to write the equations for enzyme synthesis, degradation, conversion between forms, etc. This may lead to new time dependent parameters for activating and inactivating enzymes, mRNA's etc. These new time dependent parameters can in turn be eliminated by writing the equations which govern their pool sizes. If this program is pursued ruthlessly it will eventually produce equations for all dynamically relevant chemical constituents including repressors, promoters, genetic states, 109

etc. When this process is finished all internal time dependent parameters will have been eliminated, leaving only a case of transient motion. This line of reasoning produces a criterion by which to judge models. The best model would be one with no internal time dependent parameters. The model presented is, in this sense, incomplete. Indeed, there are no practical models of biochemical differentiation which are not incomplete in this respect. For example, putting questions of experimental fact and artifact aside, timing maps of transcription and translation are nothing more than time dependent parameters and in themselves give little information about the dynamics of the process. Again it is c o m m o n to believe that all of the interesting action in differentiation occurs in the area of D N A gene state kinetics. This assumes the answer before the dynamic analysis has been performed. If a complete model were to be achieved and then simulated so as to always focus attention on the active hierarchical level, that is, the dynamic pools which are controlling the timing of the process, then this active level would move around. At some time it would display variables close to D N A metabolism, but at other times it would display variables that are far from D N A processes, such as nutrient availability during the vegetative phase of the D. discoideum life cycle. It is only when this dynamic structure is finally deciphered that we will begin to have a true understanding of differentiation. It is also commonly believed that D N A must be the major basis of control because so mucfi information is packed into DNA. But that is not where all the information is. It also takes information to specify the milieu which embeds the DNA. This may seem to be a small amount of information compared to that encoded in the DNA. However, the dynamically relevant information is less than the total information stored there. This is discussed by NEWMAN 12, where the information in use is distinguished from the information it takes to characterize a system. This was illustrated earlier where it was shown that under the steady state approximation the amount of information needed to characterize a network could be much less than under a dynamic situation. This is also shown in Figure 10. Here two genes are shown which code for two enzymes which each contain, say, 200 amino acid units. The genes then contain 110


(6O0 BASE PAIRS 1200 BITS OF I N F O R ~ T I O N )












l NO

Fig. 10. Information Shielding. Only a small portion of the genetic information is used during any given epoch. The active dynamical variables contribute a comparable amount of information.

600 base pairs and it takes 1200 bits of information to specify them. Assume that the first enzyme catalyzes the transformation of a metabolite from a pool with concentration Xl and that this reaction is always in the linear region of the enzyme kinetics. Then the reaction rate is given by klxa where kl is the rate constant for the reaction in the linear region. If it is necessary to have this rate to within 1%, then approximately 8 bits of information is needed to specify the pool concentration x I and 8 bits to specify the rate constant k 1. (This gives 1 part in 256 with some room for scale.) The 1200 bits of information at the D N A level have been reduced to 8 bits of information and this is comparable to the information needed to specify the initial state of the system. Of course, in many cases the rate law won't be so simple and it will take more information to characterize the enzymatic reaction. Yet this will still be comparable to the information it takes to specify the initial state of the system. In some cases the information will even be less. In Figure 10, gene2 codes for an enzyme which catalyzes a nonlimiting, nonbranching reaction on a pathway. In this case, only one bit of information is needed to characterize the reaction.


There is one trajectory in a high dimensional state space which describes the transformation occurring in differentiation. An adequate description of this trajectory will require the use of highly nonuniform time and space scales. The fineness of these scales will inevitably be dictated by the dynamics of the

process and is not open to free choice. The evolution of the dynamic embedding structure is an important characteristic of differentiating systems. During any given epoch there will be an active hierarchical level. This level will be composed of the dynamically active variables and will be a distinct subset of all variables in the system. This active level can be closely identified with the critical variables of the system in that they determine and limit the rate of transformation of the system during that epoch. Some techniques for simplification are as follows: a subset of variables can be studied by projecting the trajectory into subspaces. Similar subspaces can be contracted into one space producing bundles of trajectories. Metabolic networks can be reduced by using the steady state approximation and by parameterization. Time dependent parameters can be used but these represent weak points in the model and at the same time opportunities for further progress. By studying projections of individual segments of the system trajectory into various subspaces a complete description can eventually be reconstructed. This paper has been presented with all the mathematics disguised. What is really involved behind the simple descriptions is the theory of stiff differential equations and hierarchical systems, metabolic network theory and the generation of reduced equivalent networks, numerical approximation of functions, the iterative solution of equations, information theory and differential geometry. All of this and more is useful and necessary to our ultimate understanding of differentiation. The mathematics is complex - but who expected it to be simple? What is important is that the complexity be placed where it belongs, in the time evolving dynamical structure of these systems and not in a mere proliferation of variables.

References 1. Park, D. J. M., J. Theor. Biol. 46, 31-74 (1974). 2. Wright, B. E. and Gustafson, G. L., J. Biol. Chem. 247, 7875-7884 (1972). 3. Wright, B. E. and Park, D. J. M., J. Biol. Chem. (in press) (1974). 4. Rosen, R., Dynamical System Theory in Biology, Vol. 1, Wiley-Interscience, New York (1970). 5. Wright, B. E., Critical Variables in Differentiation, Prentice-Hall, Inc., Englewood Cliffs, N.J. (1973). 6. Gear, C. W., Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Inc., Englewood Cliffs, N.J. (1971). 7. Traub, J. F., Interative Methods for the Solution of Equations, Prentice-Hall, Inc., Englewood Cliffs, N.J. (1964). 8. Bastin, T., in Towards a Theoretical Biology, 2. Sketches (Waddington, C. H., ed.) Aldine Publishing Company, Chicago (1969). 9. Waddington, C. H., The Strategy of the Genes, Allen and Unwin, London (1957). 10. Goodwin, B. C., Temporal Organization in Cells, Academic Press, New York (1963). 11. Misner, C. W., Theorne, K. S. and Wheeler, J. A., Gravitation, W. H. Freeman and Company, San Francisco (1973). 12. Newman, S. A., J. Theor. Biol. 28, 411~413 (1970).

Appendix: Steady State Solution of the Flux Divider Network Pools P4(G1P) and Ps(G6P) are assumed to be in rapid equilibrium so that P5 = 6 P4- The concentration of this combined pool is given by x 1 and the concentration of P6 (UDPG) by x 2. In the rate equations for r3 and r¢ we must use P4 = xl/7 and P5 = 6xl/7. Then the equivalent network of Figure 11 is used. Reaction R 5 is suppressed but by the steady state conditions r 5 = r 3.

R~ (glycogen phosphorylase) is given by a rapid equilibrium random bi bi with B = Q model. However, this reaction is nearly saturated and the reaction rate is primarily a function of the [

P4 ~ P5



Acknowledgement This work was supported by Grants HDO4667 and HDO5357 from the National Institutes of Health, United States Public Health Service. The authors wish to acknowledge the assistance of Richard BEAUVAISin preparing the illustrations.


\ r4 = klXl/(k° ~ x2)

] [



) x2



r 3 : k2×lX 2 r 2 : k3x2

Fig. 11. The fluxdividernetworkin D. discoideurn carbohydratemetabolismused in the numericalsolution. 111

specific activity of the enzyme so rl is treated as an independent parameter. R2 (cellulose synthetase) is approximated by a linear mass action rate equation, r2 - - k 3 x 2 . R3 is actually a composite of two reaction steps, T6P synthetase and T6Pase. The second reaction is so fast that T6P has not actually been detected in D. discoideum cells. The rate of the reaction would then be given by the rate of T6P synthetase. The actual mechanism here is not known for certain. If an ordered bi bi with P, Q ~ 0 reaction mechanism is used and the measured Michaelis and inhibition constants substituted then the reaction is well approximated by a mass action rate law, r3 = k2xlx2. R4 (UDP-glucose pyrophosphorylase) is given by an ordered bi bi with P ~ 0 reaction mechanism. The reaction operates in the linear region but product inhibition by U D P G is significant so a rate equation r4 = klxl/(ko + Xz) is used. The solution to the network is obtained by setting up the steady state condition for each pool in the network and then solving for the pool sizes and reaction rates. The steady state equations are: P4 + P s : r l - klXl/(ko + X2) -- k2xlx2 = 0, P6: klXl/(ko -k x2) - k2XlX2 - k3x 2 = 0 These can be combined to give a cubic equation for xz. X2 3 -~- ClX2 2 -[- CzX 2 -~ C 3 =


where c 1 = rl/k 3 + ko, c 2 = rlko/k 3 ÷ kl/k2, c3 = - r l k l / ( k 2 k 3 ) . Once a solution for x2 is obtained then X1 = k3x2(k 0 q- x2)/(k 1 - k2x2(k 0 -+- X2)), and r2, r 3 and r4 can be immediately calculated. The cubic equation is solved by an iterative computer procedure. The best parameters available for our model are ko = 0.13 mM, kl = 1.05 mM m i n - 1, kz = 2.85 mM- 1 m i n - 1 k3 = 0.52 m i n - 1 where the concentrations are measured on a packed cell volume basis. The solution for these parameter values is given in Table 2. 112

It should be noticed that in the cubic equation for x2, and the equation for r3, the kinetic parameters kl and k2 appear only as the ratio (kl/k2). The rate r2 must also depend only o n the ratio (kl/k2) because of the steady state relation r2 = rl - 2r3. Hence an activation of R4 can be exactly offset by an equal activation of R3. It is possible to derive an expression for the limiting value of the cellulose flux, rz, as rl approaches infinity. This is given by (r2)limit


(k3/2)(-ko + x/k02 + 4(k~/k2)).

For the case calculated in Table 2 this limit is 0.284 mM/min.

Mathematical modeling of differentiation in Dictyostelium discoideum.

Methods for the dynamic analysis of biochemical differentiation are presented. These are demonstrated in the analysis of biochemical differentiation o...
1MB Sizes 0 Downloads 0 Views