The most common childhood brain cancer is medulloblastoma (MB), a devastating disease that can result in permanent serious brain damage even when successfully treated [26–28]. The etiology of MB is believed to involve aberrant activation of Sonic hedgehog (Shh) signalling [27–31]. Potential therapeutic targets are therefore regulators of the hedgehog pathway [32–34]. Although much remains to be learned about the transcriptional hierarchy involved in hedgehog activation, enough experimental data is available to construct pathway models that are quite complex. Hedgehog regulation of medulloblastomas will be used as an example to illustrate our modeling approach.
A model of hedgehog regulation of tumor growth will be constructed hierarchically, starting with a few molecular components and their role in tumor growth. Figure 1 presents the major hedgehog components that control "tumor growth". Tumor growth is a tissue-level quantity is regulated by a hierarchy of steps starting with Shh expression and hedgehog activation. Initially, we are interested in identifying the common network features or motifs that are involved in this system and building an extensible model that captures our understanding of the processes involved. The rules used to build this model are shown beside the figure. There are 8 nodes in this model and 15 reactions. Reactions for most variables include production and decay. Shh does not have production or decay reactions; it is set either to either its lowest or highest value in simulations. Specific motifs within this pathway will now be isolated and discussed for their relevance to the model. The rules for each of these cases are the same is for the whole pathway, except that lightly shaded parts are set to constant values.
Shh binds to the tumor repression gene Patched (Ptc1), releasing its suppression of Smoothened (Smo), which activates Gli transcription factors. The Gli transcription factors are the primary effectors of Shh signalling in the developing cerebellum [30, 35]. Of particular interest is regulation of the Gli family of transcription factors that are believed to be involved in medulloblastoma tumorigenesis [26, 28, 29, 33]. Several promising cancer drugs that reduce medulloblastoma growth in mouse models purportedly interfere with Gli1 and Gli2 expression by blocking the signalling pathway downstream of Smo. These include HhAntag [34] and cyclopamines [32, 36, 37]. A target of Gli1 is PTC1; Ptc1 is a transmembrane protein that blocks Smo expression when Shh binds to it. Gli1 activates Ptc1 expression and thus counters, to some extent Shh repression of Ptc1 activity [33, 35, 38]. The feedback from Gli1 to Ptc1 that is active in normal cells, as illustrated by figure 2, may play an important role in determining the effectiveness of drugs that regulate Gli1 expression downstream of Smo. Some reports suggest that this feedback loop may be inactivated in medulloblastomas in vivo, but not in cultured cells [33, 34].
Basic fuzzy network modeling concepts
We define a dynamic, interacting system by specifying the all of the variables that participate in the system and the processes (reactions) in which each node participates. The participating variables, also called nodes as in graph theory or Petri nets, can represent protein concentrations, gene expression levels, tumor size, temperature or any other variable of interest, including discrete variables.
Reactions or processes define dynamic changes to node values and must involve one or more nodes as reactants. Though the language used leans toward biochemical reactions, a useful feature of this rule-based model is the ability to connect dynamics on multiple scales. Nodes can represent any quantity for which rules of change can be defined. A reaction is defined by a reaction rate constant and a list of participating nodes together with their role in the reaction. The rate constant is a scaling factor that multiplies the reaction rate that is determined by changing reactant values. Reactions can also be defined as discrete events. Reaction rate has no meaning in this case; the event happens when the reactants satisfy the rule for firing and nothing happens when those conditions are not satisfied. This is useful for modeling processes such as chromosome separation or cytokinesis. Discrete events may also be a useful model for representing complex processes like metastasis as a single event that follows a long series of cellular reactions and sets into motion other processes.
For each reaction shown, the reactants have default rules based on their role as substrate, product, activator or inhibitor. These default rules can be overridden by giving new rules explicitly, as will be described in the methods section or may be optimized using automatic methods to integrate other data sources. The latter will not be discussed in detail here. A diagram of a system as in figure 1 is used to construct a computer model by first listing all of the participating nodes. Each reaction is then listed, with a reaction rate. Substrates are variables that are used up in a reaction; products are produced. Activators and inhibitors participate in a reaction, but their value doesn't change. The default rules for each of these roles are as follows:
When a substrate is absent, the reaction rate is zero. As the concentration of substrate increases, the reaction rate increases proportionately. Thus, if the substrate level is low, its contribution to the reaction rate is low. The reaction rate, multiplied by the rate scaling factor or rate constant, determines how fast the substrate is used up in that reaction. By default, substrates have a stoichiometry of -1. This can be changed manually as needed.
A product is produced in a reaction. It does not affect the reaction rate, except that when the product concentration reaches its highest allowable amount the reaction rate goes to zero. No more can be produced. The default stoichiometry is +1.
Activators have the same default rules as a substrate: as the concentration increases, the reaction rate contribution increases. Activator concentration is not changed by a reaction in which it participates as an activator, it only influences the rate. Activators and inhibitors have a stoichiometry of 0, indicating that the reaction does not change their concentration.
An inhibitor causes a reaction rate to slow as its concentration increases. When absent, an inhibitor has no effect on a reaction rate, which is then wholly determined by the concentrations of substrates and activators. As the amount of inhibitor increases, the reaction rate slows. The concentration at which an activator or inhibitor strongly affects a reaction rate can be set when known. By default, the effect is proportionate to the enzyme concentration.
In each reaction, the contributions of all participating reactants are averaged using one of several averaging schemes. We found that either a harmonic average or minimum rule works best. Both of these have the property that if the rule for any reactant makes the rate zero, the rate will be zero. Thus, if a reaction includes two substrates and an activator where the concentration of one substrate is zero, but the activator and the other substrate concentration are high, the output rule is zero: the reaction rate should be zero because one of the substrates is absent.
In rule-based simulations the meaning of "low" and "high" are determined by the dynamic range or universe of discourse for each variable. In the hedgehog model, most variables have a concentration range of [0.0, 1.0] for simplicity in illustrating the dynamics of particular motifs. However, the dynamic range of Shh is defined to be [0.001, 0.1] in this model. This was set to illustrate that variables can be assigned values in any desired range. Rules are defined based on the meaning of "low" and "high" within the context of the appropriate range for a given variable. Thus, a concentration of 0.1 is high for Shh in this model, while 0.1 is low for Gli1. In addition, the meaning of "low" and "high" for a given variable can be different in different reactions. Thus, in reactions where Gli1 is a product, such as a transcription reaction, low, medium and high on a scale of 0 to 1 might be represented by 0.4, 0.6 and 0.8. In a catalytic reaction where Gli1 is a potent non-stoichiometric activator, a concentration of 0.001 might be "high" in that it causes the reaction to achieve its maximum rate. This information can be put into the model rules if available. If not, model simulations can still be run using whatever information is available.
Common motifs
Positive and negative feedback
Feedback motifs are common in biological systems [39–43]. Negative feedback loops are essential for maintaining equilibrium and resisting change while positive feedback structures amplify signals to enable switch-like behaviour [39]. Gli1 is a transcription factor that promotes transcription of Ptc1 [28, 30, 33, 44]. Figure 2 highlights the Gli1 feedback loop in the hedgehog pathway. NmycP and the transforming event required for tumor growth are constitutively "on" in figure 2 simulations in order to isolate the effects of Gli1 feedback on Ptc1 and tumor growth. When Shh signalling is active, Ptc1 is repressed, causing increased Gli1 expression. When the Gli1 feedback loop is active, as in figure 2.a, the system is self-regulating and tends toward an equilibrium level of Ptc1, Smo and Gli1, as increased Ptc1 expression tends to upregulate Gli1 even in the presence of Shh [33]. When the feedback loop is broken, as in figure 2.b, Shh signalling causes rapid suppression of Ptc1 to low levels and growth of Smo and Gli1 to high levels. Disconnecting the feedback loop is accomplished in the model by simply removing Gli1 as a reactant from the production reaction for Ptc1 (rule 3 in figure 1).
Feedforward motif
Cells have a remarkable ability to orchestrate precise sequences of events using imprecise components. An important mechanism for accomplishing this feat is to filter out some signals and respond to others. A feedforward motif can filter the effects of transient signals while allowing sustained signals to activate downstream components [45]. Figure 3 isolates a feedforward motif in the hedgehog pathway. To isolate the dynamics of the feedforward motif, the transforming event was constitutively "on" and the Gli1 feedback loop to Ptc1 was deleted. Shh production was on for one-half day at the start of day 1, and then turned off. Though NmycP responded quickly to this transient signal, it was not sufficient to activate Gli1 transcription. On day 4 a longer Shh signal was turned on. Gli1 transcription began approximately one day after Gli1 production began. Tumor growth begins when Gli1 reaches a moderately high level. When Shh transcription ends, NmycP levels begin to drop immediately, but Gli1 decline is delayed. Thus, a series of prolonged aberrant Shh activation events could sustain tumor growth. Note that times and rates in this example were implemented to illustrate modeling concepts and, though physiologically reasonable, may not be accurate.
The feedforward system in figure 3 offers a potential explanation for the occurrence of medulloblastoma in haploinsufficient Ptc+/- mice. Approximately 14–20% of heterozygous Ptc+/- mice develop medulloblastoma tumors, whereas these tumors are rare in Ptc+/+ mice [33, 34]. Lower production of Ptc1 in heterozygous mice may make the probability that a perturbation to the hedgehog pathway will result in decreased Ptc1 expression for a long enough time to allow tumor growth to initiate. Transient Shh upregulation, as seen in the simulation, is not sufficient to allow sufficient Gli1 increase if Ptc1 production is high.
Single Input Motif (SIM)
In a related manner, temporal coordination of developmental processes can be achieved by differential response to a common signal [45, 46]. Single input motifs involve activation of several parallel pathways by a single activator [45]. In figure 4, Shh activates three separate pathways. Each of these parallel pathways is required for tumor growth and each has a different activation threshold, decay threshold and rate. The rules for this simulation are as for the basic model in figure 1, with the Gli1 feedback loop deleted to illustrate this motif. The simulation curves have the same shape and characteristic response profile as the differential equation simulations in [45], demonstrating that intuitive rules define the correct dynamic behaviour.
Chemical reaction kinetics
Although analytical equations are known for basic reaction kinetics and approximations for enzyme-catalyzed reactions, it is desirable to use simple rules to describe such reactions when the goal is to integrate these with other processes for which analytical equations are not known. A number of software packages are available for modeling reaction dynamics with differential equations and can be found, for example, on the systems biology workbench site [47]. Generic numerical software that might also be suitable can be found in a number of places, including the netlib numerical software library [48]. Rule-based representations of chemical reactions are relatively simple to write and easily modified, which is useful when exploring the system effects of different possible reactions through simulation.
Because of the universal approximation properties of fuzzy logic models, they are able to represent any nonlinear process as accurately as desired, given sufficient data. It may be easier to make very precise chemical reaction models with differential equations, but robust system models that integrate many levels, from approximate cellular chemistry to clinical data, may be more useful for some medical research and more easily constructed with linguistic rules by biologists and clinicians without numerical modeling experience.
Rules to determine reaction rates from substrate concentrations can be determined in a variety of ways. The mathematical details for fuzzy logic rules most appropriate for modeling chemical kinetics remain to be derived. For modeling biological systems, simplicity with reasonable accuracy is a high priority. A second order reaction and an enzyme catalyzed reaction, both involving activation of Smo by Ptc1, are implemented with fuzzy rules and examined here.
An essential part of hedgehog pathway activation involves release of Ptc1 repression of Smo. Smo activates downstream transcription factors of the Gli family, which in turn upregulate cell cycle components [31, 33, 38, 49]. One model suggests that Smo is in a balance between active and inactive states. Smo is activated when it combines stoichiometrically with a small molecule (sm) that changes its conformation. The small molecule must be actively transported across the cell membrane by Ptc1, which acts catalytically.
Figure 5 shows a model of Smo reaction with a small molecule agonist in a second order reaction after transport facilitated by Ptc1. Analytical solutions to first, second and third order kinetics are compared to simulations from a rule-based model for Smo activation by a small molecule. The simulated time series for Smo activation follows second order kinetics quite closely. Analytic solutions for first, second and third order kinetics are shown for comparison. This is intended to show that relatively simple rule-based models give reasonable biochemical dynamics. We note here that the default rules for first order reactions, where the rate is proportional to substrate concentration only, is mathematically equivalent to the usual differential equation for a first order reaction. Third order reactions are also well-approximated by simple reaction rules with three substrates.
The model in figure 5 is simulated with the following rules, with only three essential reactions. For each reaction, take the harmonic mean of substrate concentrations (zero/v. low/low/med/high/v. high); then reaction rate is (zero/v. low/low/med/high/v. high). The three reactions are: Reaction 1: Smo_inactive and sm are substrates and active Smo is the product. The rate coefficient is 1.0e-3 s-1. Reaction 2: the reverse reaction, has one substrate, activated Smo-sm complex, and two products, with the rate coefficient 1.0e-3 s-1. Reaction 3: The sm is transported across the membrane at a rate much faster than it is used by Smo and is continually supplied at a constant rate that is catalyzed by Ptc1.
To analyze the system dynamics of the hedgehog pathway, it may be sufficient to model the catalytic nature of Ptc1 action heuristically using rules that capture the nonlinear effects of Ptc1 on Smo: if Ptc1 is zero, Smo activation is at its normal or highest level, Vmax. When a small amount of Ptc1 is present, Smo concentration is only 20% of its highest level (a linguistic value of "very low"). If Ptc1 is higher than that, the equilibrium concentration of active Smo goes to near zero. Normally, the linguistic terms zero, very low, low, medium, high and very high are assumed to mean fractions 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, respectively, of the maximum possible value. These are used as defaults for all rules. However, these can be adjusted when more information is available. In figure 6, Ptc1 catalyzes Smo inactivation. It was reported that a concentration of Ptc1 of 1/45 the concentration of Smo reduces the level of active Smo by 80% [50]. The meaning of very low, low and medium we set at 0.012, 0.022 (1/45) and 0.032 for Ptc1 when it acts enzymatically to deactivate Smo. The rule for this reaction was: the inactivation rate is zero when Ptc1 is zero or very low and very high when Ptc1 is low or higher. The curve for equilibrium Smo versus Ptc1 concentration in figure 6.b shows the typical S-shape for an enzymatic reaction, with the greatest change in rate occurring for Ptc1 concentration of "low", the value at which the rate goes from zero to very high.
It is important to emphasize that the goal here is not just to show that fuzzy rules can be found to match differential equation models of chemical kinetics. Rather, a model based on intuitively reasonable rules that describe the reaction qualitatively produce reaction dynamics that are very close to dynamics derived from chemical kinetics using rigorous mathematical derivations.
Oscillators, switches and discrete events in the cell cycle
Cell cycle machinery plays a central role in cell proliferation, apoptosis and cell fate determination. Failure to exit the cell cycle at the correct time may cause abnormal development and aberrant re-entry into the cell cycle by terminally differentiated cells may be a primary cause of many cancers. Understanding regulation of the cell cycle is therefore a central issue in the application of molecular biology to medicine. For that reason, a number of differential equation models of the cell cycle in yeast have been constructed for research purposes [51–54]. Regulation of the mammalian cell cycle is more complex than that of yeast, though many homologous genes and proteins and common network motifs are involved.
A model of the mammalian cell cycle was constructed with several important components. Chromosome separation and cell division were included as discrete events that occur when certain conditions are met. The cell cycle is an oscillatory system that involves not only continuous processes, but a series of discrete events [55]. Discrete events are an important mechanism in cells for controlling the precise timing of key events and insuring irreversibility with imprecise components [56, 57]. Thus, even though many cell cycle parameters are poorly constrained, the overall behaviour of the cell cycle system is well characterized and heuristic modeling approaches are suitable for studying the control dynamics or testing possible ideas for new drug targets [58].
Figure 7 presents a cell cycle model that illustrates how chemical kinetics, enzyme catalyzed reactions and discrete events, such as chromosome separation, cell division, and entry/exit to a combined G1SG2 phase can be integrated. The primary molecular components of this model are cyclinD1, cyclinB2, APC, and Cdc14. Cell cycle models with these or similar components have been constructed using differential equations [54, 59].
Tumor growth is linked to this cell cycle via cell division. Tumor size in this case is a measure of the number of cell divisions that have occurred. If the cell cycle arrests, tumor size shrinks slowly through a continuous apoptotic reaction. Rules for this cell cycle model are as follows:
1. cyclinD1 production is continuous
2. cyclinD1 decay is proportional to cyclinD1 and APC levels and is rapid
3. cyclinB production is proportional to cyclinD1 level
4. cyclinB decay is zero when APC is less than high, is low when APC is high, and is very high when APC is very high.
5. Cdc14 production is proportional to cyclinB levels
6. Cdc14 decay is continuous and proportional to Cdc14 concentration
7. APC production is proportional to Cdc14 level, but is zero when Cdc14 is very low and low.
8. APC decay is directly proportional to cyclinB, cell mass and entry into G1SG2 and determined by the harmonic average of these.
9. Chromosome separation event occurs when cyclinD1 is medium or higher and G1SG2 is in the on state (G1SG2 is a discrete event, either on or off). G1SG2 turns off when chromosome separation occurs.
10. Cell division occurs when cyclinB is low or less and decreasing at any rate; chromosome separation event must also have occurred. G1SG2 turns on when cell division occurs and chromosome separation turns off.
11. Cell division turns off when cell mass is very low or less and cell division is currently on.
12. Cell mass growth is constant. When cell division occurs (a single event), cell mass is halved. Cell mass has a range of values between one and two.
The growth rate of cell mass has a controlling effect on the rest of the cycle. It was set to make the cycle length approximately 1 day, a typical value for human cells undergoing mitosis [60]. Figure 8 shows time courses for molecular components and events in the cell cycle model, as well as tumor growth. The shape and ordering of curves for cyclinB, Cdc14 and APC are very similar here to those shown in Tyson's cell cycle model [59]. Chromosome separation is a discrete event that is on for a brief time at the beginning of mitosis and then turns off when cell division occurs. Another discrete event, cell division, is not shown here. If either of these events does not occur (turn on), the cell cycle halts. Similarly, knocking out any of the critical components of this pathway, cyclinD1, cyclinB, Cdc14 or APC, for example by deleting their production reaction, will also cause the cycle to halt. In figure 8.b the cell cycle was set to stop suddenly on day 20. When this happens, tumor growth, a measure of the number of cell divisions that have occurred, stops and the tumor slowly decays through the apoptotic reaction. The latter is a simple reaction with tumor size as the sole participant (as a substrate that is used up).
Modeling and drug target discovery
Computational models are potentially very useful for exploring the effects of manipulating complex pathways through drugs or other means. Rule based models of two features that make them promising tools for drug target discovery. First, it is easy for experts to use their domain knowledge to build models. Rules can be changed easily without mathematical manipulations to carry out thought experiments with many hundreds or thousands of components. Secondly, algorithms for automatically manipulating fuzzy rules have been developed for a variety of engineering applications and include genetic algorithms and neural nets. It may be possible for experts to build first-stage models of pathways and then apply computational methods to optimize the models using other data sources, or to search for optimal drug targets.
To illustrate the use of the model for carrying out computational experiments a rule-based model was constructed using elements of the previously presented pathways. Our goal in this example is to explore the effects of a drug, HhAntag, that reportedly interferes with Gli1 activation of cyclinD1, thereby stopping progression of the cell cycle [34]. HhAntag was introduced as a new component to the model in figure 1. HhAntag inhibits both Gli1 and, to a lesser extent, Gli2 expression downstream of Smo.
An interesting reported effect of HhAntag on cultured neural tissues is that a relatively low dose suppresses Gli1 levels, but a thousand fold increase is needed to halt tumor growth [34]. This effect is not observed in vivo. An important key to this difference is the suggestion that the feedback loop from Gli1 to Ptc1 is operates in the cultured system but not in vivo. As discussed previously, Ptc1 acts catalytically on Smo expression: very low levels of Ptc1 are sufficient to significantly reduce Smo activity [29]. As shown in figure 9, when the Gli1 feedback loop to Ptc1 is active small doses of drug are sufficient to suppress tumor growth (figure 9.b). Without this feedback, the drug is not as effective in low doses. Only at much higher doses can the drug suppress both Gli proteins sufficiently to stop tumor growth (figure 9.c).
It should be emphasized here that the hedgehog pathway is not well understood and the model presented here is in preliminary stages. Nevertheless, it illustrates important modeling ideas. The next step is to replace the tumor growth module with a cell cycle model. First, we develop and test a preliminary cell cycle model using rules. This will then be connected to the hedgehog pathway.
A more complete model for hedgehog regulation of the cell cycle is shown in figure 10. The cell cycle model discussed previously has been combined with the components in figure 9. Integration of these rule-based models was simple and required little more than pasting the rule files for each into one file and modifying the reactions for cyclinD1 so that Gli1 and Gli2 regulated its production. In this model, Shh directly induces Nmyc expression and indirectly affects Nmyc posttranslational modification, mediated by its indirect targets, cyclinB and possibly other cyclins [30, 44, 61]. Shh also binds to Ptc1, as described previously, activating Gli expression. The tumor growth process from figure 8 is now replaced by the cell cycle model from of figure 9. Note that a circular pathway exists among activated cyclinB, Nmyc and its phosphorylated form NmycP, and cyclinD1. An additional component or transforming event, mentioned in figure 1, is needed to start this cycle, which reinforces the suggestion that another yet unknown target may also be required for cell transformation [61].
Once started, the cycle will continue as long as the Shh pathway remains active. Mitotic degradation of Nmyc permits neuronal precursor cell cycle exit in the absence of Shh signalling or in the case of an intrinsic program-directed shift toward differentiation. Thus neural fate specification and cancer activation are regulated by the same machinery.
In figure 10.b simulation time courses show the effect of decreasing insulin growth factor (IGF) on the cell cycle. In this figure, IGF is decreased beginning on day 6, allowing GSK-3β to increase. Nmyc-P is then phosphorylated and degraded. Note that Nmyc-P initially decreases, and then recovers for one cell cycle through complicated feedbacks. Incorporating as many of these complicated interactions into a model is essential for identifying potential drug targets, as cells are remarkably robust and have built-in fault tolerance systems that may not be evident from examination of static network or interaction diagrams [56, 62, 63]. The effects of drug dose will be the same for this model as in figure 9, as the primary action of the drug on Gli2 has not been changed. However, adding details to the cell cycle now allows more detailed investigation of the interaction between the drug, its targets, and cell cycle components.