A cellular automaton model of cellular signal ... - Semantic Scholar
Computers in Biology and Medicine 30 (2000) 1±21
www.elsevier.com/locate/compbiomed
A cellular automaton model of cellular signal transduction Jens U. Wurthner a,*, 1, Amal K. Mukhopadhyay a, Claus-JuÈrgen Peimann b a
Institute for Hormone and Fertility Research, University of Hamburg, Grandweg 64, D-22529, Hamburg, Germany b Institute for Computerscience in Medicine, University of Hamburg, Martinistr. 52, D-20246, Hamburg, Germany Received 18 December 1997; accepted 17 September 1999
Abstract On the basis of cellular automata models, a software speci®cally tailored to model biochemical reactions involved in cellular signal transduction was implemented on a personal computer. Recent data regarding desensitization processes in mouse Leydig cells are used to simulate the underlying reactions of signal transduction. Pretreatment of real Leydig cells with dierent molecules results in a modi®cation of the signal transduction cascade leading to a diminished response of the cells during subsequent stimulations. The model is capable of simulating the complex behavior of this intracellular second messenger production in a qualitative and semi-quantitative way. The results indicate that quantitative simulations on a molecular level will be possible once appropriate computer hardware is available. The simulations and results of the cellular automaton presented are easily described and comprehended, which make it a useful tool that will facilitate research in cellular signal transduction and other ®elds covering complex reaction networks. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: Cellular automata; Desensitization; Mouse Leydig cells; Signal transduction; Simulation
1. Introduction The standard-technique used to model natural phenomena is to set up mathematical equations describing the system (mostly partial dierential equations). Often these equations can only be approximated by means of numerical integration through computers. This * Corresponding author. Tel: +1-301-594-5729; fax: +1-301-596-8395. E-mail address: [email protected] (J.U. Wurthner). 1 Present address: Laboratory of Cell Regulation and Carcinogenesis, National Cancer Institute, Bldg. 41 Rm. C629, 41 Library Drive, Bethesda, MD 20982, USA. 0010-4825/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 0 - 4 8 2 5 ( 9 9 ) 0 0 0 2 0 - 7
2
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
approach is rather tedious, complicated and time-consuming. In addition, in many cases a substantial computational power is needed. Furthermore, the results of the computations can be dicult to interpret. An alternative approach is to utilize Cellular Automata (CA). Their value for simulations of natural phenomena has been discussed previously (for a comprehensive review, see [1,2]). Fundamental aspects of nature, locality and parallelism leading to complex behavior, are immanent features of CA. Also, hardware optimized for CA runs exceedingly fast and oers computational power at a fraction of the cost if compared to normal computers [3]. CA have been successfully employed for many years in various ®elds. Applications in the biomedical sciences were scarce until recently, when they became useful to study complex biological systems, for example the immune system [4,5], arti®cial, self-replicating systems [6], tumor growth [7] and excitation in cardiac tissue [8]. We have now devised a CA for the simulation of biochemical networks, especially signal transduction pathways. The biochemical data generated in this ®eld are already too enormous and diverse for single researchers to comprehend. Computer programs that can simulate those reaction networks could be a valuable tool to help integrate the information about cellular signaling pathways, to discover interconnections between pathways and to help generate hypotheses and design meaningful experiments to prove them. To test whether CAs are capable of modeling signal transduction pathways, we decided to simulate transmembrane signal transduction in mouse Leydig cells, which is a well characterized experimental system. In short, activation of a transmembrane signaling process is initiated by the hormones gonadotrophin luteotropin (LH) or chorionic gonadotropin (CG) following their interaction with cell surface receptors. This leads to an increase of the intracellular concentration of the second messenger cyclic adenosine monophosphate (cAMP) and a subsequent stimulation of steroid-hormone biosynthesis [9]. The stimulation of steroidhormone biosynthesis can be reduced by stimulating the cells repeatedly with either the two hormones mentioned above [10,14], or with other molecules that act through dierent pathways, such as phorbol esters like 4-b-phorbol 12-myristate 13-acetate (PMA) [11]. Another molecule, phosphodiesterase (PDE), is an intracellular enzyme that plays a role in ending a cAMP response by hydrolyzing cAMP molecules. We now show that the signaling pathway as well as the more sophisticated mechanisms of desensitization can be simulated by our model and lead to the generation of results closely resembling biochemical data. The reader should be aware that studies performed in the laboratory will be addressed as experiments and experiments done with the model as simulations in the following text. Also, the term cells, which usually describes pixels in the CA literature, is not used in this way here. Cells in this manuscript refer to biologic cells that are used in experiments or simulated with the CA.
2. Background Most cellular automata, including the one presented here, use pixels of various colors (also called states) on a computer screen in a simulation window as representatives for objects that are to be simulated. Therefore, in the context of cellular automata, the term pixel does not
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
3
merely refer to a dot on a computer screen, but to a real object that is represented by that pixel. Which object is represented by such a pixel is speci®ed in a table which assigns each pixel color to a speci®c object. For example, green pixels can be assigned to sulfuric acid molecules, blue pixels to water molecules etc. Pixels can be viewed either in their entity, in which they form a complex pattern in the simulation window, or in equal sized subunits, in which they form so called neighborhoods. In CA, a neighborhood is a spatial arrangement of pixels in which each participating pixel has an in¯uence on one or more pixels within that spatial arrangement. Fig. 1 depicts the neighborhood used in the model presented in this paper, the Margolus-Neighborhood [3,12]: by de®nition, the spatial arrangement is a two by two array of pixels. Prior to the simulation, each possible pixel combination of the two by two pixel array Ð in this paper referred to as source block Ð is assigned to a user-de®ned combination of four pixels, referred to as the destination block. In the simplest case, depicted in Fig. 2, there are two pixel states, white and black, leading to 16 dierent combinations (source blocks) that are assigned to 16 combinations speci®ed by the user (destination blocks). The assignment of source blocks to destination blocks is a necessary prerequisite for a simulation. Since the states are internally represented by numbers, the list of source and destination blocks is in fact stored as an array of numbers, as exempli®ed by Table 1. The state of each pixel in the array is coded by a digit in a speci®c position in the entire number that represents the combination. The whole table of source blocks and assigned destination blocks is called a rule table, as it contains the rules that govern the replacement of pixels and, therefore, the simulation. It is stored in the random access memory. In the case of the example above, the 16 dierent destination blocks can be manually entered by the user. In the case of more elaborate simulations with larger rule tables, as in our model, special algorithms have to be implemented that ®ll in the rule table based on more general descriptions of reactions, e.g. chemical reactions. This eliminates the need to individually ®ll in each entry of a large rule table. The number of individual entries in this rule table, i.e. the number of possible pixel combinations, depends on the number of pixels in a neighborhood and the number of states each pixel can
Fig. 1. Margolus Neighborhood. The Margolus Neighborhood is by de®nition an array of four pixels arranged in a square (two by two grid). During a simulation, the screen is divided into a multitude of such squares, which are replaced one by one according to speci®c rules stored in a rule table.
4
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
Fig. 2. Diusion rule for Margolus Neighborhood. A graphical representation of a rule that can be utilized by a cellular automaton based on the Margolus Neighborhood to simulate diusion processes. The blocks on the left side of each arrow constitute so called source blocks, the ones on the right side destination blocks. During a simulation, each source block from the screen is compared to the source blocks in the rule table, and there will be exactly one match. The assigned destination block is then taken by the algorithm to replace the source block on the screen. A rule table contains all possible combinations of pixels on the side of the source blocks, so for every source block found on the screen there is exactly one entry in the rule table. In this case, two dierent pixel colors lead to 16 dierent pixel combinations.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
5
Table 1 Representation of the rule table in memory. Each individual reaction is stored in a coded format in consecutive bytes of computer memory. Each of these reactions is represented by four bits describing the source block and four bits describing the destination block. Table 1 corresponds to the graphical representation of the same rule that is displayed in Fig. 2 Combination number
Source block
Destination block
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0000 I000 0I00 II00 00I0 I0I0 0II0 III0 000I I00I 0I0I II0I 00II I0II 0III IIII
0000 000I 00I0 00II 0I00 0I0I I00I 0III I000 0II0 I0I0 I0II II00 II0I III0 IIII
have. For the Margolus-Neighborhood, therefore, the number of possible combinations equals (number of pixel states)4, which results in 160000 combinations in the case of the 20 dierent states we chose for our model. Increasing the size of the neighborhood will not only exponentially increase the memory demand but will also increase the computation time, so that a two by two array serves as a good compromise. As implied by the name, source blocks are replaced by destination blocks during a simulation. This is achieved by the reaction algorithm which scans the simulation window and sequentially picks up source blocks, i.e. four pixel arrays, searches for the same source block in the rule table, looks at the destination block assigned to this source block, and then replaces the source block with the destination block in the simulation window. Theoretically, all the two by two arrays in the simulation window can be replaced simultaneously, since they are independent of each other, which explains why parallel computing is ideally suited for CA. With the PC used for the model described here, however, the replacement of pixel arrays was done in a sequential fashion. In the CA described herein, the user can decide if a molecule is stationary or moving. Moving molecules always move in a diagonal fashion for the following reason: after each two by two array in the simulation window is replaced, one arti®cial time unit called generation is said to be passed. During the following generation, the virtual grid that overlies the simulation window and determines the upper left corner of each two by two array is shifted one pixel to the right and one pixel downward. Every pixel that had been in the lower right corner of a two by two array is now in the upper left corner, every pixel that had been in the upper right
6
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
corner in the lower left corner etc. In the next generation, the grid is shifted back to its original place. This grid shifting has one important eect; it allows motion and thereby enables diusion and reaction of molecules that were not within the same neighborhood at the beginning of the simulation. An example depicted in Fig. 3 clari®es this concept: a pixel in an upper left corner of a two by two array is replaced by a destination block that contains this pixel in the lower right corner (reaction number 2 of Fig. 2). After this replacement, the grid moves, which puts this pixel of a lower right corner again in an upper left corner. This will again lead to a replacement according to the same rule, putting the pixel in a lower right corner. Then the grid shifts back, and the pixel is again in an upper left corner. If no other pixels interfere, the pixel will move across the whole simulation window in a diagonal fashion, from the upper left corner towards the lower right corner. Depending on their initial placement, single pixels will therefore move in one of the four possible diagonal directions. The diagonal fashion of movement is speci®ed in reactions 2, 3, 5 and 9 of the rule in Fig. 2. With a two by two grid, only four directions of movement are possible. Movements of straight up, down, left and right are also possible but not used with the CA discussed here. In the case of the graphical representation of a rule table, depicted in Fig. 2, the simulation leads to a global behavior of diusion of the white and black pixels [1]. The interpretation of what white and black pixels stand for is irrelevant for the simulation process as such. However, the rules entered by the user usually re¯ect how the pixel states are to be interpreted. In addition, we
Fig. 3. Virtual grid enables pixel motility. After each source block has been replaced in a generation (1), the virtual grid that determined where the source blocks are in relation to the pixels on the screen, moves one pixel downward and one pixel to the right (2). Then the next generation of replacement follows (3). Without introduction of a virtual grid that moves back and forth every generation, the rules lined out in Fig. 2 would lead to an oscillation of pixels, but not to movement and diusion. For example, a pixel in the upper left corner would be replaced according to rule number 2 in Fig. 2 by a block with a pixel in the right lower corner. Without grid movement, this block would be subject to rule number 9 in Fig. 2 in the following generation, which again places the pixel in the upper left corner of the same two by two array, which is its former position. However, with grid movement Ð as shown in the middle diagram Ð the pixel is positioned in an upper left corner of a source block, which is replaced again according to rule number 2 in Fig. 2. Numbering of the gridlines is intended to facilitate the perception of pixel movement. The solid arrow indicates the direct replacement of the block, the hatched arrow indicates the passing of a whole generation, i.e. the replacement of all other blocks on the screen, before the grid moves.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
7
chose to provide an input mask in which the user enters the meaning of the pixel states, i.e. the molecule names. The result of a simulation is generally represented by either the global arrangement of pixels in the simulation window and its changes over time Ð as in the example of a simple diusion process Ð or by the change in the number of pixels with a certain state, e.g. representing a concentration change of a certain molecule. The latter form of result is utilized in our model generating Figs. 5, 8 and 9. 3. Material and methods Based on the Margolus neighborhood [3,12] shown in Fig. 1, a software was programmed that extended the capabilities of the original algorithm in the following ways: It allows for up to 20 dierent molecules to be simulated, it allows for molecules to be both motile and non motile, it allows input of molecules and reactions in a format used in biology or chemistry without abstraction, generating the more abstract rules for the rule table automatically, it incorporates the use of reaction constants, and it automatically generates results in the form of molecule counts (which represent molecule concentrations). Before a simulation can be performed, four requirements have to be met. 1. 2. 3. 4.
The The The The
user has to enter the names of the molecules. user has to enter the reactions. user has to enter the starting pattern of molecules on the screen. program has to generate a rule table.
3.1. The molecules In a molecule table that holds 20 entries, the user speci®es which molecules are to be simulated, which of those are to be moving during a simulation, and which ones are not. For example, cell membrane molecules are normally immotile (otherwise cells would disintegrate), whereas molecules such as hormones are denoted ``motile'' so that they can diuse to dierent regions of the simulation window. By entering the molecules into the table, each molecule is automatically assigned a unique state. Pixels on the screen of that state represent the corresponding molecule in a simulation. Table 2 depicts such a molecule table. The assignment of molecule names to states facilitates the input of reactions, since it allows the use of actual chemical equations instead of states. For the cellular automaton it is irrelevant what names are assigned to what states, i.e. how the user interprets the pixels in the simulation window. Therefore, this information is not used when the rule table is generated, and consequently it does not in¯uence the simulation as such. 3.2. The reactions In a reaction table that can hold several thousand entries, the user speci®es which of his
8
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
Table 2 Molecule table. The user assigns a molecule to each color that can be used by the cellular automaton and whether pixels of this color are motile or not. The motile property is used during rule generation, and the color information during the actual simulation when the simulation is displayed in a simulation window Number
Molecule
Mobility
Represented by color
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
background cell membrane human choriogonadotrophin inactive receptor activated receptor inactive G protein active G protein inactive adenylate cyclase active adenylate cyclase cyclic AMP phosphodiesterase phosphodiesterase-inhibitor inhibited phosphodiesterase phosphodiesterase-cAMP-complex inactive protein kinase C active protein kinase C desensitized G protein phorbolester AlFÿ 4 not de®ned
yes no yes no no no no no no yes yes yes yes yes no no no yes yes yes
white dark yellow purple dark blue blue black grey orange green red light blue magenda dark purple pink light green light orange dark gray cyan light yellow light pink
previously entered molecules react with which, and what the resulting molecules are of such reactions. Forward and backward reactions are considered distinct and have to be entered as separate reactions. The user also determines the likelihood with which reactions occur, i.e. the reaction constant. This reaction constant is proportional to the association constant if the reaction describes an association, or is proportional to the dissociation constant if the reaction describes a dissociation. For example: A reaction B+R 4 BR has an association constant of 1000 l/mole. This has to be arbitrarily converted into a reaction constant of e.g. 15000, as there is no conversion rule of how to transform an association constant to a reaction constant of the CA per se. Another reaction with an association constant of e.g. 2000 l/mole now has to be assigned a value twice that assigned for the reaction with 1000 l/mole: 30000. This demonstrates that the ®rst association constant in a given set of simulated reactions has to be assigned an arbitrary reaction constant, and the other association constants have to be assigned reaction constants that re¯ect their ratio to the ®rst constant. The range 0 to 32000 of possible values for the reaction constant is a compromise between memory-requirement (15 bit) and accuracy. A range that covers natural reaction constants would have been more desirable but was not possible due to the same memory constrictions limiting the number of dierent molecules.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
9
3.3. The starting pattern of molecules on the screen A typical simulation window contains one or more cells, which is/are represented by a square composed of cell membrane pixels. The interior of the cell is initially ®lled with background and intracellular molecules, e.g. phosphodiesterase. The cell membranes are two pixels in diameter to allow for the incorporation of signal transducing units composed of receptors, G-proteins and adenylate cyclase proteins. Constructing the cells and the complexes by hand, i.e. by drawing with the mouse in the simulation window, is dicult and prone to error. Therefore, special automatic drawing algorithms are implemented. The composition and number of the signal transduction units to be placed into the membrane is user de®ned, and in the simulations described here it contains an inactive receptor, an inactive G-protein and an inactive adenylate cyclase (AC) at the beginning of a simulation, and sometimes also an inactive protein kinase C molecule. Fig. 4a and b depict two examples. As recent evidence obtained by our group [13] points to the possibility that both CG and PMA induced desensitization is mediated via a common, protein kinase C dependent pathway, the model integrates an activation of protein kinase C as part of the desensitization process. Although it is thought that protein kinase C is translocated from the cytosol to the plasma membrane during this desensitization, it is directly incorporated into the signal transduction units for reasons of simplicity (Fig. 4b). However, initially the protein kinase C is inactive, and only after activation it can promote desensitization. 3.4. The generation of a rule table The rule table is based on the Margolus neighborhood and contains 10 bytes per rule. The ®rst 4 bytes contain the source block, 1 byte per pixel, the second 4 bytes contain the destination block, also 1 byte per pixel, and bytes 9 and 10 contain the reaction constant for that particular rule. Out of the ®rst and second four bytes, each byte no. 1 represents the pixel in the upper left corner, each byte no. 2 the pixel in the upper right corner, each byte no. 3 the pixel in the lower right corner, and each byte no. 4 the pixel in the lower left corner. Per default, the whole rule table is ®lled with rules that will lead to diagonal movement of all molecules when the software is initialized. This is achieved by calculating all possible combinations and ®lling them into the source blocks of the rule table. The destination blocks are subsequently calculated by rotating bytes 1 through 4 of the source block two times. This will put byte no. 1 (upper left pixel) into byte no. 3 (lower right pixel), no. 2 into no. 4 etc. Empirical data (not shown) demonstrated that in the case of only two molecules that diagonally interact, a rotation by only one byte instead of two leads to a faster distribution of molecules on the simulation window. Therefore, the calculation of default destination blocks was altered accordingly. In CA that allow only two possible states, 16 possible rules exist, which are depicted in Fig. 2 (graphical representation) and Table 1 (coded form). In the case of the CA described herein, 20 dierent molecules lead to 160000 possible rules. Normally, the reactions entered by the user in a non coded form constitute only a subset of those 160000 rules. For example, the reaction ``A+B+C+D 4 E+F+G+H'' is encoded by 24 dierent rules, since there are exactly 4 3 2 1 possibilities for the molecules ``A'', ``B'', ``C'' and ``D'' to interact on a two by two grid. The program calculates these permutations and the
10
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
corresponding permutations of the product molecules ``E'', ``F'', ``G'' and ``H''. With ``A'' being represented by 1, ``B'' by 2 and so forth, the ®rst rule of this reaction is therefore ``1,2,3,4:5,6,7,8'', the second one ``1,2,4,3:5,6,8,7'', and so forth, all re¯ecting dierent possible spacial arrangements of four dierent pixels in a two by two grid. It is crucial to note that the ®rst four bytes of these rules were calculated before at initialization time and used to build the rule table. By generating the rules from the reactions, the program only has to calculate the dierent permutations of the molecules of the reactions entered by the user and then ®nd these
Fig. 4. (a) Signal transduction unit. A signal transduction unit is de®ned by the user, each individual molecule is speci®ed by assigning a molecule name from the molecule list entered previously to each molecule in the unit. The relevant molecules for the signal transduction in this unit are the inactive receptor, the inactive G-protein and the inactive adenylate cyclase (AC). This information is used by an algorithm that automatically draws cells on the screen: Depending on the number of signal transduction units wanted per cell, the algorithm calculates the number of cell membrane pixels needed in between the units and draws the cells accordingly (the cell membrane molecules will be extended to the left and the right of each unit in order to generate a cell with a continous border). The cell membrane molecules within the transduction unit are necessary to avoid holes in the membrane at the signal transduction sites. (b) Signal transduction unit including protein kinase C. This unit is similar to the one depicted in Fig. 4a, except that it contains an additional molecule, inactive protein kinase C (PKC). The fact that the inactive receptor is incorporated in the membrane instead of protruding does not aect the signaling. In order to make the inactive PKC accessible for modifying molecules (such as phorbolesters) both from outside and from within the cell, a background ``molecule'' was incorporated into the membrane. This allows for intracellular molecules to diuse right next to the PKC, which is necessary for an interaction.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
11
permutations in the rule table. In the above example, these are ``1,2,3,4'', ``1,2,4,3'', and so forth. The program will ®nd these 24 permutations in the rule table and then write the corresponding permutations of the product molecules ``E'', ``F'', ``G'' and ``H'', i.e. ``5,6,7,8'', ``5,6,8,7'', etc., into the 4 bytes reserved for the destination block, overwriting the default values. When the CA encounters the more common form of two molecules reacting with each other, e.g. a hormone-molecule attaching to its receptor, the program calculates the permutations by adding two other molecules to complete the four byte code (i.e. two by two array). These unspeci®ed molecules can be any of the available molecules. For example, assuming that hormone is encoded by ``3'' and receptor by ``5'', the left side of the equation ``hormone+receptor 4 active receptor+background'' will be represented by ``3, 5, 1, 1'', ``3, 5, 1, 2'', ``3, 5, 1, 3'', . . . ``3, 5, 19, 20'' and ``3, 5, 20, 20'' and sought in the rule table (the total number of combinations in this case is 4 3 202=4800). The corresponding destination blocks will then be overwritten with 4 bytes coding for one activated receptor state (coded as 6), one background state (coded as 0), and two other states corresponding to the ones calculated for the source block, e.g. ``6,0,1,1'', ``6,0,1,2'' . . . The reason for the right side of the equation to contain a distinct ``background'' is that this is a de®ned state which should not be replaced with any other molecule by the rule generating algorithm. After these four requirements are met, the user can start the simulation. An algorithm now scans the simulation window and picks up every two-by-two pixel source block in that window, arranges those 4 pixels in a row of 4 bytes, and then looks for the same 4 bytes (=source block) in the rule table. The corresponding 4 bytes of the destination block are subsequently used to replace the two by two pixels in the simulation window. The simulation algorithm also determines a random number between 0 and 32000 for every source block it encounters. Only if the reaction constant previously entered by the user is greater or equal to this number, the replacement of the source block by the destination block will take place. Scanning and replacing each block in the whole simulation window in that fashion is thought to take place simultaneously (even though it is a sequential process on non-parallel computers), and the time span of one update to the next represents the smallest unit of time in a CA, one ``generation''. All molecules, regardless of mass and volume, are represented by states of pixels in the simulation window. Moving molecules like cAMP and hormones travel in a diagonal fashion, with a change of direction after collisions with other molecules. Molecules that are thought to be not moving as much in relation to moving ones, such as membrane molecules, are ®xed in a certain position that is determined initially by the user de®ned starting pattern. The model was implemented on an Apple Macintosh Computer running system 7.x. Hyper Card (Claris, Cupertino, CA, USA) and Think C version 5.0 (Symantec, Cupertino, CA, USA) were used as programming languages.
4. Results Initial simulations demonstrated that the model comprises essential aspects of bimolecular reactions, such as hyperbolic concentration kinetics, proportionality between reaction constants and reaction kinetics, the observation of the law of mass eect and ®rst and second order
12
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
kinetics (data not shown). Another important characteristic of biologic systems, the presence of concentration gradients and their diminution over time, could be modeled by placing a multitude of molecules in a particular area of the simulation window. During the subsequent simulation, the molecules moved around the simulation window, and since there were fewer molecules to re¯ect from in the direction of the concentration gradient, the molecules on average traveled in that direction, thereby leveling the concentration gradient. If desired, membrane molecules could be de®ned in a way that promote only the passage of some molecules but not others, simulating a semipermeable membrane. Impermeable membranes Ð as in the case of the following simulations Ð are simulated by membranes that re¯ect all moving molecules. In the simulation represented by Fig. 5, a transmembranous signal transduction process is simulated. It reproduces a typical cAMP kinetic of mouse Leydig cells stimulated with CG, which is shown in Fig. 6 for comparison. As outlined in the Methods Section, the following four requirements were met before the simulation was started: The user entered the molecules, shown in Table 2, the reactions, shown in Table 3, and the starting pattern, which showed six cells with a total of 304 inactive signal transduction complexes automatically incorporated into
Fig. 5. Simulation of intracellular cAMP kinetics. Two cells with a total of 304 inactive signal-transduction units as shown in Fig. 4a are stimulated with 100 molecules of CG. Intracellular cAMP degradation is carried out by 100 phosphodiesterase (PDE) molecules, which are in turn partly inhibited by 300 molecules of PDE inhibitor. The underlying reactions as entered by the user are shown in Table 3, the rule table automatically generated based on these reactions contains 160000 individual rules (data not shown).
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
13
Fig. 6. Gonadotrophin induced intracellular cAMP-kinetics in `real' mouse-Leydig cells. Mouse Leydig cells (50000 per 0.5 ml) were incubated with two dierent concentrations of hCG (R 5 ng/ml, * 20 ng/ml) in the presence of a phosphodiestarese inhibitor (0.25 mM IBMX) for dierent periods of time at 378C. The cells were then separated from the medium by short centrifugation (1 min at 10000 g), and the determination of intracellular cAMP was carried out with a speci®c radioimmunoassay. Results are from three independent determinations; vertical bars represent 2SD, which is omitted where it is smaller than the symbol.
the cell membranes, identical to the one shown in Fig. 4a. There were also 100 CG molecules placed into the simulation window outside the cell boundaries, and 100 phosphodiesterase molecules and 300 phosphodiesterase inhibitor molecules placed within the cell boundaries. The CA then calculated the rule table with its 160000 individual rules. The simulation was then done in a window of 240 320 pixels for 750 generations, which constitutes a compromise between speed and accuracy of simulation and is therefore used in most of the simulations. The essential characteristics of the experimentally derived kinetics is captured by the simulated kinetics, such as an initially high rise in cAMP concentration, a peak value of cAMP (near to generation 350), and a slow degradation of the cAMP. This slow degradation was accomplished by introducing a phosphodiesterase and its inbitor into the model. The phosphodiesterase degrades cAMP, but it is slowed by the inhibitor. The graph shown in Fig. 5 was obtained by automatically counting the pixels with the state that represents cAMP molecules every 10 generations and then by plotting these counts against the generations. As the generations are a form of arti®cial time, such a plot resembles a time course of the concentration of whatever molecule was counted.
14
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
Table 3 Reactions that constitute the signal transduction rule used to generate the results of Fig. 5. On the left side the educts are listed, on the right side the products; reaction constants are 32000 except for reactions 2 (c = 4000) and 11 (c = 1000) Educt side of the reactions 1. hCGa+inactive receptor 2. background+active receptor 3. active receptor+inactive G-protein 4. active G-protein+inactive ACc 5. inactive receptor+active G-protein 6. active G-protein+active AC 7. active AC+cellmembrane+background+background 8. cAMPb+phosphodiesterase 9. cAMP-phosphodiesterase-complex 10. phosphodiesterase+inhibitor 11. inhibited phosphodiesterase
Product side of the reactions 4 4 4 4 4 4 4 4 4 4 4
background+active receptor hCG+inactive receptor active receptor+active G-protein active G-protein+active AC inactive receptor+inactive G-protein inactive G-protein+active AC active AC+cellmembrane+cAMP+cAMP cAMP-phosphodiesterasecomplex+background phosphodiesterase inhibited phosphodiesterase phosphodiesterase+inhibitor
a
CG: human choriogonadotrophin. cAMP: cyclic adenosine monophosphate. c AC: adenylate cyclase.
b
Fig. 7 demonstrates how the ®rst steps of signal transduction are actually simulated by the model. It is important to recognize that only a very small part of the whole simulation window is shown on an enlarged scale. Initially, all transducing molecules within the membrane are inactive. When a hormone diuses into the vicinity of a receptor, this receptor becomes activated. The hormone molecule disappears, because it is thought to be part of the active receptor (reaction no. 1 in Table 3). The active receptor now activates the G-protein, which ®nally activates the adenylate cyclase. Depending on the reaction constants, the adenylate cyclase molecules become inactive again after one or more generations. When the receptor becomes inactive, the hormone appears again and diuses away (generation x + 2 in Fig. 7). In order to incorporate ampli®cation, which is one of the principals of signal transduction, one or more of the transducing molecules will revert more slowly to the inactive state compared to how fast it became activated. Also, each active adenylate cyclase Ð at least in the rules underlying the simulation presented in this paper Ð produces two cAMP molecules per generation (generation x + 4 in Fig. 7, reaction no. 7 in Table 3). The ®nal simulations presented here focus on homologous and heterologous desensitization. This desensitization is achieved by simulating the pre-incubation of six cells either with the hormone CG (homologous, Fig. 8) or with the phorbol ester PMA (heterologous, Fig. 9). According to the reactions entered (similar to those shown in Table 3), CG activates the receptor, which activates the G-Protein, which activates adenylate cyclase as well as protein kinase C. In contrast, PMA directly activates protein kinase C. Activated protein kinase C in turn modi®es the G-protein so that it can not be activated by an active receptor (`desensitized G-protein'). The pre-incubation with CG spans about 1200 generations, the one with PMA about 600 generations. In both cases about 180 out of 240 G-proteins are desensitized after the pre-incubations, the dierent time spans being a result of dierent reactions constants involved
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
15
Fig. 7. (Caption overleaf).
16
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
in the activation of protein kinase C. Prior to re-stimulation, the simulation windows are cleared of all protein kinase C, CG and PMA molecules (simply by deleting the corresponding pixels). Re-stimulations are undertaken with either CG or AlFÿ 4 . The results depicted in Fig. 8 show clearly that pre-incubation with CG leads to a diminished response in terms of cAMPaccumulation during a second stimulation with either CG or AlFÿ 4 . A similar result is observed after pre-incubation with PMA (Fig. 9). The standard deviations were obtained by recording cAMP kinetics in three dierent simulations that contained the same cells and 100 molecules of CG. For each simulation, the CG molecules were distributed in a dierent random manner on the simulation window. 5. Discussion In this study we have described a computer model for cellular signal transduction pathways based on a cellular automaton. Although CA have been used before to model natural phenomena, they were mainly applied in the ®eld of physics, i.e. ¯uid dynamics. In recent years there were some attempts to utilize the advantages of CA in the medical ®eld, focusing on immunology, cardiology and experimental oncology. However, models to simulate signal transduction cascades have had a mathematically based approach so far, using dierential equations and special algorithms to solve them [15±17]. To our knowledge, this is the ®rst report that characterizes a CA speci®cally designed for the purpose of cellular signal transduction. The ®rst step in building a model is to clearly formulate a theory of how the dierent parts of the system to be studied interact. The CA model presented here requires the input of the molecules which play a part in the signal transduction, their concentrations, and the reaction constants of all reactions involved. Not all of this information is presently known, so some assumptions had to be made. For example, the stoichiometry of receptors to G-proteins and to adenylate cyclase molecules was arbitrary chosen to be 1:1:1 in the signal transduction units. Arbitrary numbers were also taken for other molecules, such as hormone molecules, since it was the purpose to show qualitative behavior of the cellular automaton rather than exact quantitative results. Arbitrary in this context means the absolute number of molecules of each kind taking part in the simulation, not the number of molecules that interact in a reaction, which is ®xed in the rules as exempli®ed in Table 3. The reaction constants that were used were based on our own experimental observations and
Fig. 7. Simulation of signal transduction. This sequel depicts a small part of the simulation window in the course of 5 generations during a simulation of signal transduction. A small part of the cellmembrane that contains a single signal transduction unit as explained in Fig. 5a is shown. Once a CG pixel is in direct contact with the receptor, the sourceblock that contains the receptor and the CG molecule is substituted with a destination block that contains only an activated receptor, leaving a background pixel where the CG had been. In the following generation (x + 2), this leads to activation of the G-protein, and at the same time to inactivation of the receptor, which releases the CG. In generation x + 3, the active G-proteins turns inactive again, and the adenylate cyclase becomes active. This leads to the production of two pixels of cyclic Adenosine Monophosphate (cAMP), and the inactivation of the adenylate cyclase in generation x + 4.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
17
Fig. 8. Simulation of adenylate cyclase activity after CG induced desensitization. Six cells with a total of 240 signaltransduction units as displayed in Fig. 2b were incubated for 1200 generations with (.) or without (w) 100 molecules of CG in a 184 320 pixels large simulation window. Afterwards all CG and PKC molecules were removed, and a restimulation was carried out with 100 molecules of CG (a) or AlFÿ 4 (b). The number of cAMP molecules measured are means of three dierent simulations 2SD.
published data of signaling cascades. Since the possible numbers 0 to 32000 cannot fully span all reaction constants possible in such systems, their values were assigned on a rational basis with regard to known constants and then modi®ed by trial and error. However, only a minimal amount of trial and error was necessary, which involved only the correction of reaction constants in order to have the model generate results in a reasonable time. The fact that there were no corrections to be made on the algorithms of the CA nor the fundamental reactions of
Fig. 9. Simulation of adenylate cyclase activity after PMA induced desensitization. The simulation setup was as described for Fig. 5, except that cells were ``preincubated'' with 100 molecules of phorbolester prior to removal of this PKC-stimulating agent. Restimulation was carried with CG (a) or AlFÿ 4 (b).
18
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
the model suggest that this is due to the naturalistic way CA emulate nature, i.e. generation of complex behavior based on simple rules and their application to a large number of interacting molecules. Our extension of the CA model by incorporating random numbers that in¯uence the basic reactions takes deterministic properties out of these low level calculations, whereby it most likely enhances the generation of complex behavior. A number of simpli®cations of the system were necessary to formulate the model. This was partly due to limitations of a simulation as such, and partly due to the limitations forced upon the CA by the computer hardware. Simpli®cations made in our CA include the principal restriction of all structures like cells to two dimensions, the restricted possibility of molecule movement in only four dierent directions in these two dimensions, and the restriction of only four dierent molecules interacting at one time in a given neighborhood (because it is restricted to a two by two array). In addition, characteristics such as molecular weight, molecular speed and other physical properties cannot be simulated. Furthermore, the total number of molecules being simulated is limited to 475000 at the moment, which leads to problems when very small numbers of molecules are supposed to react with each other, for example a reaction of ATP with an activated adenylate cyclase. If the simulation were to be carried out in the way that an activated adenylate cyclase would only form a cAMP molecule if it had contact with an ATP molecule, the whole simulation window would have to be ®lled with ATP to guarantee a reasonable production of cAMP. Obviously this is not possible, so simpli®cations had to be introduced to avoid these kinds of problems. In addition, due to its limitations in the number of possible reaction constants and the arbitrary selection of the ®rst reaction constant in a given set of reactions, a simple relation of the backward reaction constants and the forward reaction constants with the equilibrium constant does not exist. Nevertheless, it was demonstrated that other basic laws of chemical reactions, such as exponential reaction kinetics, the law of mass eect and the law of mass conservation are obeyed by the model. Despite these limitations, simulations were fairly accurate in the behavior that was modeled, even when only a small number of hormone molecules were used (100 molecules of CG in the simulation depicted in Figs. 8 and 9). Repeated simulations with a dierent spatial distribution of those molecules resulted in similar cAMP kinetics, as evidenced by small standard deviations of cAMP concentration at various time points. One might expect with only a few CG molecules stochastic errors of small samples, resulting in dierent cAMP kinetics. Apparently, the numbers of molecules used here were already large enough to minimize these errors. Based on experimental observations, our model ful®lls a number of important general characteristics of cellular signal transduction: . . . . . .
the hormone binds speci®cally to its receptor the hormone±receptor-complex activates a G-protein the active G-protein activates other eector molecules, such as adenylate cyclase the eector molecule produces a second messenger (in this case cAMP) there is an ampli®cation of the signal, expressed as a low hormone/cAMP ratio there is an activation of both second messenger generation and desensitization
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
19
The simulation results of the CA presented here show that CA are capable of dealing with simulations of biochemical reaction cascades, especially those involved in cellular signal transduction, in a qualitative and semi-quantitative way (semi-quantitative in the sense that reaction constants can be used to generate dierent reaction-kinetics in a network of reactions). In order to approach simulations on a true molecular scale that will predict exact concentrations of molecules at dierent time points, CA have to be further developed, mainly with respect to size and speed. Parallel to this development, knowledge about the components involved in the cellular signaling cascades will increase, enabling us to incorporate more detail into the models. Since the number of reactions and components is limited, ultimately all of these will be known. At that point the bene®t of simulations will become more obvious, because the pure knowledge of the reaction mechanisms will not be sucient to predict and thoroughly understand cellular behavior. The data presented in this paper show that CA based on the methods which we developed can serve as a basis for the development of systems that will be able to perform such simulations.
6. Summary On the basis of cellular automata (CA) models, a software speci®cally tailored to model biochemical reactions involved in cellular signal transduction was implemented on a personal computer. Although in recent years there were some attempts to utilize the advantages of CA in the biomedical ®eld, models that simulated signal transduction cascades had only utilized dierential equations and special algorithms to solve them. To our knowledge, this is the ®rst report that characterizes a cellular automaton speci®cally designed for the purpose of cellular signal transduction. Initial simulations with this model demonstrated that it comprises essential aspects of bimolecular reactions, such as hyperbolic concentration kinetics, proportionality between reaction constants and reaction kinetics, the observation of the law of mass eect and ®rst and second order kinetics. We then tailored the model for the simulation of transmembrane signal transduction in mouse Leydig cells, a process which is known to involve the hormones gondadotrophin luteotropin and/or chorionic gonadotropin. Following their interaction with cell-surface receptors, these hormones are capable of activating a transmembrane signaling process that leads to an increase of the intracellular concentration of the second messenger cyclic adenosine monophosphate (cAMP) and a subsequent stimulation of steroid-hormone biosynthesis. However, pretreatment of Leydig-cells with dierent stimulators results in a modi®cation of the signal transduction cascade leading to a diminished response of the cells during subsequent stimulations, a process known as desensitization. Recent data regarding these desensitization processes in Leydig cells were subsequently used to simulate the underlying reactions of signal transduction. The model is capable of simulating this complex behavior of the intracellular second messenger production, which involves the activation of second messenger generation and desensitization processes. Moreover, it ful®lls other important biological characteristics of cellular signal transduction, namely speci®city of hormone±receptor interactions, the interaction of signaling intermediates, and the ampli®cation of the extracellular signal, expressed as a hormone/intracellular-messenger ratio. The simulation results of the cellular automaton presented here show that CA are capable of
20
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
simulating biochemical reaction cascades, especially those involved in cellular signal transduction, in a qualitative and semi-quantitative way and suggest that quantitative simulations on a molecular level will be possible once appropriate computer hardware is available. As simulations and results of the CA are easily described and comprehended, CA are useful tools that will facilitate research in cellular signal transduction and other ®elds covering complex reaction networks. Acknowledgements The authors would like to thank Prof. Leidenberger for moral and ®nancial support to carry out this work, and Jan WuÈrthner for many programming tips. References [1] T. Tooli, N. Margolus, Cellular Automata Machines Ð A New Environment for Modeling, MIT Press, Massachusetts, 1990. [2] H. Gutowitz, Cellular Automata Ð Theory and Experiment, MIT Press, Cambridge, Massachusetts, 1991. [3] N. Margolus, T. Tooli, Cellular automata machines, in: E. Doolen (Ed.), Lattice Gas Methods for Partial Dierential Equations, SFI SISOC, Addison-Wesley, Cambridge, Massachusetts, 1990, pp. 219±249. [4] H.B.A. Sieburg, Logical dynamic systems approach to the regulation of antigen-driven lymphocyte stimulation, in: H.B.A. Sieburg (Ed.), Theoretical Immunology, Addison-Wesley, 1988, pp. 273±291. [5] H.B. Sieburg, J.A. McCutchan, O.K. Clay, L. Cabalerro, J.J. Ostlund, Simulations of HIV infection in arti®cial immune systems, Physica D 45 (1990) 208±227. [6] J.A. Reggia, S.L. Armentrout, H. Chou, Y. Peng, Simple systems that exhibit self-directed replication, Science 259 (1993) 1282±1287. [7] A.S. Qi, X. Zheng, C.Y. Du, B.S.A. An, Cellular automaton model of cancerous growth, J. Theor. Biol. 161 (1993) 1±12. [8] D.W. Frazier, P.D. Wolf, J.M. Wharton, A.S. Tang, W.M. Smith, R.E. Ideker, Stimulus induced critical point. Mechanism for electrical initiation of reentry in normal canine myocardium, Journal of Clinical Investigation 83 (1989) 1039±1052. [9] J.M. Saez, Leydig cells: endocrine, paracrine, and autocrine regulation, Endocr. Rev. 15 (1994) 574±626. [10] M. Schumacher, M. Schwarz, W. BraÈndle, Desensitation of the cAMP system in mouse Leydig cells by hCG, cholera toxin, dibutyryl cAMP: Localization of the `lesion' to the guanine nucleotide regulatory protein± adenylate cyclase complex, Mol. Cell. Endocrinol. 34 (1984) 67±80. [11] A.K. Mukhopadhyay, M. Schumacher, Inhibition of hCG-stimulated adenylate cyclase in puri®ed mouse Leydig cells by the phorbol ester PMA, FEBS Lett. 187 (1985) 56±60. [12] T. Tooli, N. Margolus, Programmable matter: Concepts and realization, Physica D 47 (1991) 263±272. [13] J.U. WuÈrthner, M. Kistler, M. Kratzmeier, A.K. Mukhopadhyay, LH/hCG-receptor is coupled to both adenylate cyclase and protein kinase C signaling pathways in isolated mouse Leydig cells, Endocrine J. 3 (1995) 579±584. [14] M. Schumacher, G. SchaÈfer, V. Lichtenberg, H. Hilz, Maximal steroidogeneic capacity of mouse leydig cells, FEBS Lett. 107 (1979) 398±402. [15] J.E.A. McIntosh, R.P. McIntosh, Mathematical Modelling and Computers in Endocrinology, Springer-Verlag, New York, 1980. [16] M.D. Houslay, The use of selective inhibtitors and computer modelling to evaluate the role of speci®c high anity cyclic AMP phosphodiesterases in the hormonal regulation of hepatocyte intracellular cyclic AMP concentratrions, Cellular Signaling 2 (1989) 85±98. [17] S. Swillens, D. Mercan, Computer simulation of a cytosolic calcium oscillator, Biochem. J. 271 (1990) 835±838.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
21
Jens U. Wurthner obtained his MD from the University of Hamburg Medical School in 1995. In the same year he began his post graduate training in internal medicine at the University of Freiburg Medical Center in the department for hematology/oncology. In 1997 he received a DFG scholarship for cancer research and subsequently joined the Laboratory of Cell Regulation and Carcinogenesis at the National Cancer Institute in Bethesda, where he is currently working on intermediates and modi®ers of the TGF-b signaling cascade. Amal Mukhopadhyay obtained his Ph.D. in Biochemistry from the Banaras H. University, Varanasi, India in 1971 and joined this University as a Lecturer of Physiology in the same year. In 1982 he moved to Hamburg to join the Institute for Hormone and Fertility Research at the University of Hamburg, Hamburg, Germany where is currently working as the Director, Division of Reproductive Endocrinology. His area of specialization is molecular basis of hormone-induced signalling mechanism in reproductive tissues. Claus J. Peimann received his diploma in physics from the University of Hamburg, Germany in 1971 and his Ph.D. in Mathematics in 1974, also University of Hamburg. He obtained his Habilitation in Medical Informatics in 1986 at the University of Hamburg Medical School and is currently an assistant professor at the Institute for Computerscience in Medicine, University of Hamburg Medical School. His special interests include simulation methods like Petri Nets and cellular automata, Hypertext/ Hypermedia and teachware.
www.elsevier.com/locate/compbiomed
A cellular automaton model of cellular signal transduction Jens U. Wurthner a,*, 1, Amal K. Mukhopadhyay a, Claus-JuÈrgen Peimann b a
Institute for Hormone and Fertility Research, University of Hamburg, Grandweg 64, D-22529, Hamburg, Germany b Institute for Computerscience in Medicine, University of Hamburg, Martinistr. 52, D-20246, Hamburg, Germany Received 18 December 1997; accepted 17 September 1999
Abstract On the basis of cellular automata models, a software speci®cally tailored to model biochemical reactions involved in cellular signal transduction was implemented on a personal computer. Recent data regarding desensitization processes in mouse Leydig cells are used to simulate the underlying reactions of signal transduction. Pretreatment of real Leydig cells with dierent molecules results in a modi®cation of the signal transduction cascade leading to a diminished response of the cells during subsequent stimulations. The model is capable of simulating the complex behavior of this intracellular second messenger production in a qualitative and semi-quantitative way. The results indicate that quantitative simulations on a molecular level will be possible once appropriate computer hardware is available. The simulations and results of the cellular automaton presented are easily described and comprehended, which make it a useful tool that will facilitate research in cellular signal transduction and other ®elds covering complex reaction networks. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: Cellular automata; Desensitization; Mouse Leydig cells; Signal transduction; Simulation
1. Introduction The standard-technique used to model natural phenomena is to set up mathematical equations describing the system (mostly partial dierential equations). Often these equations can only be approximated by means of numerical integration through computers. This * Corresponding author. Tel: +1-301-594-5729; fax: +1-301-596-8395. E-mail address: [email protected] (J.U. Wurthner). 1 Present address: Laboratory of Cell Regulation and Carcinogenesis, National Cancer Institute, Bldg. 41 Rm. C629, 41 Library Drive, Bethesda, MD 20982, USA. 0010-4825/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 0 - 4 8 2 5 ( 9 9 ) 0 0 0 2 0 - 7
2
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
approach is rather tedious, complicated and time-consuming. In addition, in many cases a substantial computational power is needed. Furthermore, the results of the computations can be dicult to interpret. An alternative approach is to utilize Cellular Automata (CA). Their value for simulations of natural phenomena has been discussed previously (for a comprehensive review, see [1,2]). Fundamental aspects of nature, locality and parallelism leading to complex behavior, are immanent features of CA. Also, hardware optimized for CA runs exceedingly fast and oers computational power at a fraction of the cost if compared to normal computers [3]. CA have been successfully employed for many years in various ®elds. Applications in the biomedical sciences were scarce until recently, when they became useful to study complex biological systems, for example the immune system [4,5], arti®cial, self-replicating systems [6], tumor growth [7] and excitation in cardiac tissue [8]. We have now devised a CA for the simulation of biochemical networks, especially signal transduction pathways. The biochemical data generated in this ®eld are already too enormous and diverse for single researchers to comprehend. Computer programs that can simulate those reaction networks could be a valuable tool to help integrate the information about cellular signaling pathways, to discover interconnections between pathways and to help generate hypotheses and design meaningful experiments to prove them. To test whether CAs are capable of modeling signal transduction pathways, we decided to simulate transmembrane signal transduction in mouse Leydig cells, which is a well characterized experimental system. In short, activation of a transmembrane signaling process is initiated by the hormones gonadotrophin luteotropin (LH) or chorionic gonadotropin (CG) following their interaction with cell surface receptors. This leads to an increase of the intracellular concentration of the second messenger cyclic adenosine monophosphate (cAMP) and a subsequent stimulation of steroid-hormone biosynthesis [9]. The stimulation of steroidhormone biosynthesis can be reduced by stimulating the cells repeatedly with either the two hormones mentioned above [10,14], or with other molecules that act through dierent pathways, such as phorbol esters like 4-b-phorbol 12-myristate 13-acetate (PMA) [11]. Another molecule, phosphodiesterase (PDE), is an intracellular enzyme that plays a role in ending a cAMP response by hydrolyzing cAMP molecules. We now show that the signaling pathway as well as the more sophisticated mechanisms of desensitization can be simulated by our model and lead to the generation of results closely resembling biochemical data. The reader should be aware that studies performed in the laboratory will be addressed as experiments and experiments done with the model as simulations in the following text. Also, the term cells, which usually describes pixels in the CA literature, is not used in this way here. Cells in this manuscript refer to biologic cells that are used in experiments or simulated with the CA.
2. Background Most cellular automata, including the one presented here, use pixels of various colors (also called states) on a computer screen in a simulation window as representatives for objects that are to be simulated. Therefore, in the context of cellular automata, the term pixel does not
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
3
merely refer to a dot on a computer screen, but to a real object that is represented by that pixel. Which object is represented by such a pixel is speci®ed in a table which assigns each pixel color to a speci®c object. For example, green pixels can be assigned to sulfuric acid molecules, blue pixels to water molecules etc. Pixels can be viewed either in their entity, in which they form a complex pattern in the simulation window, or in equal sized subunits, in which they form so called neighborhoods. In CA, a neighborhood is a spatial arrangement of pixels in which each participating pixel has an in¯uence on one or more pixels within that spatial arrangement. Fig. 1 depicts the neighborhood used in the model presented in this paper, the Margolus-Neighborhood [3,12]: by de®nition, the spatial arrangement is a two by two array of pixels. Prior to the simulation, each possible pixel combination of the two by two pixel array Ð in this paper referred to as source block Ð is assigned to a user-de®ned combination of four pixels, referred to as the destination block. In the simplest case, depicted in Fig. 2, there are two pixel states, white and black, leading to 16 dierent combinations (source blocks) that are assigned to 16 combinations speci®ed by the user (destination blocks). The assignment of source blocks to destination blocks is a necessary prerequisite for a simulation. Since the states are internally represented by numbers, the list of source and destination blocks is in fact stored as an array of numbers, as exempli®ed by Table 1. The state of each pixel in the array is coded by a digit in a speci®c position in the entire number that represents the combination. The whole table of source blocks and assigned destination blocks is called a rule table, as it contains the rules that govern the replacement of pixels and, therefore, the simulation. It is stored in the random access memory. In the case of the example above, the 16 dierent destination blocks can be manually entered by the user. In the case of more elaborate simulations with larger rule tables, as in our model, special algorithms have to be implemented that ®ll in the rule table based on more general descriptions of reactions, e.g. chemical reactions. This eliminates the need to individually ®ll in each entry of a large rule table. The number of individual entries in this rule table, i.e. the number of possible pixel combinations, depends on the number of pixels in a neighborhood and the number of states each pixel can
Fig. 1. Margolus Neighborhood. The Margolus Neighborhood is by de®nition an array of four pixels arranged in a square (two by two grid). During a simulation, the screen is divided into a multitude of such squares, which are replaced one by one according to speci®c rules stored in a rule table.
4
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
Fig. 2. Diusion rule for Margolus Neighborhood. A graphical representation of a rule that can be utilized by a cellular automaton based on the Margolus Neighborhood to simulate diusion processes. The blocks on the left side of each arrow constitute so called source blocks, the ones on the right side destination blocks. During a simulation, each source block from the screen is compared to the source blocks in the rule table, and there will be exactly one match. The assigned destination block is then taken by the algorithm to replace the source block on the screen. A rule table contains all possible combinations of pixels on the side of the source blocks, so for every source block found on the screen there is exactly one entry in the rule table. In this case, two dierent pixel colors lead to 16 dierent pixel combinations.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
5
Table 1 Representation of the rule table in memory. Each individual reaction is stored in a coded format in consecutive bytes of computer memory. Each of these reactions is represented by four bits describing the source block and four bits describing the destination block. Table 1 corresponds to the graphical representation of the same rule that is displayed in Fig. 2 Combination number
Source block
Destination block
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0000 I000 0I00 II00 00I0 I0I0 0II0 III0 000I I00I 0I0I II0I 00II I0II 0III IIII
0000 000I 00I0 00II 0I00 0I0I I00I 0III I000 0II0 I0I0 I0II II00 II0I III0 IIII
have. For the Margolus-Neighborhood, therefore, the number of possible combinations equals (number of pixel states)4, which results in 160000 combinations in the case of the 20 dierent states we chose for our model. Increasing the size of the neighborhood will not only exponentially increase the memory demand but will also increase the computation time, so that a two by two array serves as a good compromise. As implied by the name, source blocks are replaced by destination blocks during a simulation. This is achieved by the reaction algorithm which scans the simulation window and sequentially picks up source blocks, i.e. four pixel arrays, searches for the same source block in the rule table, looks at the destination block assigned to this source block, and then replaces the source block with the destination block in the simulation window. Theoretically, all the two by two arrays in the simulation window can be replaced simultaneously, since they are independent of each other, which explains why parallel computing is ideally suited for CA. With the PC used for the model described here, however, the replacement of pixel arrays was done in a sequential fashion. In the CA described herein, the user can decide if a molecule is stationary or moving. Moving molecules always move in a diagonal fashion for the following reason: after each two by two array in the simulation window is replaced, one arti®cial time unit called generation is said to be passed. During the following generation, the virtual grid that overlies the simulation window and determines the upper left corner of each two by two array is shifted one pixel to the right and one pixel downward. Every pixel that had been in the lower right corner of a two by two array is now in the upper left corner, every pixel that had been in the upper right
6
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
corner in the lower left corner etc. In the next generation, the grid is shifted back to its original place. This grid shifting has one important eect; it allows motion and thereby enables diusion and reaction of molecules that were not within the same neighborhood at the beginning of the simulation. An example depicted in Fig. 3 clari®es this concept: a pixel in an upper left corner of a two by two array is replaced by a destination block that contains this pixel in the lower right corner (reaction number 2 of Fig. 2). After this replacement, the grid moves, which puts this pixel of a lower right corner again in an upper left corner. This will again lead to a replacement according to the same rule, putting the pixel in a lower right corner. Then the grid shifts back, and the pixel is again in an upper left corner. If no other pixels interfere, the pixel will move across the whole simulation window in a diagonal fashion, from the upper left corner towards the lower right corner. Depending on their initial placement, single pixels will therefore move in one of the four possible diagonal directions. The diagonal fashion of movement is speci®ed in reactions 2, 3, 5 and 9 of the rule in Fig. 2. With a two by two grid, only four directions of movement are possible. Movements of straight up, down, left and right are also possible but not used with the CA discussed here. In the case of the graphical representation of a rule table, depicted in Fig. 2, the simulation leads to a global behavior of diusion of the white and black pixels [1]. The interpretation of what white and black pixels stand for is irrelevant for the simulation process as such. However, the rules entered by the user usually re¯ect how the pixel states are to be interpreted. In addition, we
Fig. 3. Virtual grid enables pixel motility. After each source block has been replaced in a generation (1), the virtual grid that determined where the source blocks are in relation to the pixels on the screen, moves one pixel downward and one pixel to the right (2). Then the next generation of replacement follows (3). Without introduction of a virtual grid that moves back and forth every generation, the rules lined out in Fig. 2 would lead to an oscillation of pixels, but not to movement and diusion. For example, a pixel in the upper left corner would be replaced according to rule number 2 in Fig. 2 by a block with a pixel in the right lower corner. Without grid movement, this block would be subject to rule number 9 in Fig. 2 in the following generation, which again places the pixel in the upper left corner of the same two by two array, which is its former position. However, with grid movement Ð as shown in the middle diagram Ð the pixel is positioned in an upper left corner of a source block, which is replaced again according to rule number 2 in Fig. 2. Numbering of the gridlines is intended to facilitate the perception of pixel movement. The solid arrow indicates the direct replacement of the block, the hatched arrow indicates the passing of a whole generation, i.e. the replacement of all other blocks on the screen, before the grid moves.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
7
chose to provide an input mask in which the user enters the meaning of the pixel states, i.e. the molecule names. The result of a simulation is generally represented by either the global arrangement of pixels in the simulation window and its changes over time Ð as in the example of a simple diusion process Ð or by the change in the number of pixels with a certain state, e.g. representing a concentration change of a certain molecule. The latter form of result is utilized in our model generating Figs. 5, 8 and 9. 3. Material and methods Based on the Margolus neighborhood [3,12] shown in Fig. 1, a software was programmed that extended the capabilities of the original algorithm in the following ways: It allows for up to 20 dierent molecules to be simulated, it allows for molecules to be both motile and non motile, it allows input of molecules and reactions in a format used in biology or chemistry without abstraction, generating the more abstract rules for the rule table automatically, it incorporates the use of reaction constants, and it automatically generates results in the form of molecule counts (which represent molecule concentrations). Before a simulation can be performed, four requirements have to be met. 1. 2. 3. 4.
The The The The
user has to enter the names of the molecules. user has to enter the reactions. user has to enter the starting pattern of molecules on the screen. program has to generate a rule table.
3.1. The molecules In a molecule table that holds 20 entries, the user speci®es which molecules are to be simulated, which of those are to be moving during a simulation, and which ones are not. For example, cell membrane molecules are normally immotile (otherwise cells would disintegrate), whereas molecules such as hormones are denoted ``motile'' so that they can diuse to dierent regions of the simulation window. By entering the molecules into the table, each molecule is automatically assigned a unique state. Pixels on the screen of that state represent the corresponding molecule in a simulation. Table 2 depicts such a molecule table. The assignment of molecule names to states facilitates the input of reactions, since it allows the use of actual chemical equations instead of states. For the cellular automaton it is irrelevant what names are assigned to what states, i.e. how the user interprets the pixels in the simulation window. Therefore, this information is not used when the rule table is generated, and consequently it does not in¯uence the simulation as such. 3.2. The reactions In a reaction table that can hold several thousand entries, the user speci®es which of his
8
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
Table 2 Molecule table. The user assigns a molecule to each color that can be used by the cellular automaton and whether pixels of this color are motile or not. The motile property is used during rule generation, and the color information during the actual simulation when the simulation is displayed in a simulation window Number
Molecule
Mobility
Represented by color
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
background cell membrane human choriogonadotrophin inactive receptor activated receptor inactive G protein active G protein inactive adenylate cyclase active adenylate cyclase cyclic AMP phosphodiesterase phosphodiesterase-inhibitor inhibited phosphodiesterase phosphodiesterase-cAMP-complex inactive protein kinase C active protein kinase C desensitized G protein phorbolester AlFÿ 4 not de®ned
yes no yes no no no no no no yes yes yes yes yes no no no yes yes yes
white dark yellow purple dark blue blue black grey orange green red light blue magenda dark purple pink light green light orange dark gray cyan light yellow light pink
previously entered molecules react with which, and what the resulting molecules are of such reactions. Forward and backward reactions are considered distinct and have to be entered as separate reactions. The user also determines the likelihood with which reactions occur, i.e. the reaction constant. This reaction constant is proportional to the association constant if the reaction describes an association, or is proportional to the dissociation constant if the reaction describes a dissociation. For example: A reaction B+R 4 BR has an association constant of 1000 l/mole. This has to be arbitrarily converted into a reaction constant of e.g. 15000, as there is no conversion rule of how to transform an association constant to a reaction constant of the CA per se. Another reaction with an association constant of e.g. 2000 l/mole now has to be assigned a value twice that assigned for the reaction with 1000 l/mole: 30000. This demonstrates that the ®rst association constant in a given set of simulated reactions has to be assigned an arbitrary reaction constant, and the other association constants have to be assigned reaction constants that re¯ect their ratio to the ®rst constant. The range 0 to 32000 of possible values for the reaction constant is a compromise between memory-requirement (15 bit) and accuracy. A range that covers natural reaction constants would have been more desirable but was not possible due to the same memory constrictions limiting the number of dierent molecules.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
9
3.3. The starting pattern of molecules on the screen A typical simulation window contains one or more cells, which is/are represented by a square composed of cell membrane pixels. The interior of the cell is initially ®lled with background and intracellular molecules, e.g. phosphodiesterase. The cell membranes are two pixels in diameter to allow for the incorporation of signal transducing units composed of receptors, G-proteins and adenylate cyclase proteins. Constructing the cells and the complexes by hand, i.e. by drawing with the mouse in the simulation window, is dicult and prone to error. Therefore, special automatic drawing algorithms are implemented. The composition and number of the signal transduction units to be placed into the membrane is user de®ned, and in the simulations described here it contains an inactive receptor, an inactive G-protein and an inactive adenylate cyclase (AC) at the beginning of a simulation, and sometimes also an inactive protein kinase C molecule. Fig. 4a and b depict two examples. As recent evidence obtained by our group [13] points to the possibility that both CG and PMA induced desensitization is mediated via a common, protein kinase C dependent pathway, the model integrates an activation of protein kinase C as part of the desensitization process. Although it is thought that protein kinase C is translocated from the cytosol to the plasma membrane during this desensitization, it is directly incorporated into the signal transduction units for reasons of simplicity (Fig. 4b). However, initially the protein kinase C is inactive, and only after activation it can promote desensitization. 3.4. The generation of a rule table The rule table is based on the Margolus neighborhood and contains 10 bytes per rule. The ®rst 4 bytes contain the source block, 1 byte per pixel, the second 4 bytes contain the destination block, also 1 byte per pixel, and bytes 9 and 10 contain the reaction constant for that particular rule. Out of the ®rst and second four bytes, each byte no. 1 represents the pixel in the upper left corner, each byte no. 2 the pixel in the upper right corner, each byte no. 3 the pixel in the lower right corner, and each byte no. 4 the pixel in the lower left corner. Per default, the whole rule table is ®lled with rules that will lead to diagonal movement of all molecules when the software is initialized. This is achieved by calculating all possible combinations and ®lling them into the source blocks of the rule table. The destination blocks are subsequently calculated by rotating bytes 1 through 4 of the source block two times. This will put byte no. 1 (upper left pixel) into byte no. 3 (lower right pixel), no. 2 into no. 4 etc. Empirical data (not shown) demonstrated that in the case of only two molecules that diagonally interact, a rotation by only one byte instead of two leads to a faster distribution of molecules on the simulation window. Therefore, the calculation of default destination blocks was altered accordingly. In CA that allow only two possible states, 16 possible rules exist, which are depicted in Fig. 2 (graphical representation) and Table 1 (coded form). In the case of the CA described herein, 20 dierent molecules lead to 160000 possible rules. Normally, the reactions entered by the user in a non coded form constitute only a subset of those 160000 rules. For example, the reaction ``A+B+C+D 4 E+F+G+H'' is encoded by 24 dierent rules, since there are exactly 4 3 2 1 possibilities for the molecules ``A'', ``B'', ``C'' and ``D'' to interact on a two by two grid. The program calculates these permutations and the
10
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
corresponding permutations of the product molecules ``E'', ``F'', ``G'' and ``H''. With ``A'' being represented by 1, ``B'' by 2 and so forth, the ®rst rule of this reaction is therefore ``1,2,3,4:5,6,7,8'', the second one ``1,2,4,3:5,6,8,7'', and so forth, all re¯ecting dierent possible spacial arrangements of four dierent pixels in a two by two grid. It is crucial to note that the ®rst four bytes of these rules were calculated before at initialization time and used to build the rule table. By generating the rules from the reactions, the program only has to calculate the dierent permutations of the molecules of the reactions entered by the user and then ®nd these
Fig. 4. (a) Signal transduction unit. A signal transduction unit is de®ned by the user, each individual molecule is speci®ed by assigning a molecule name from the molecule list entered previously to each molecule in the unit. The relevant molecules for the signal transduction in this unit are the inactive receptor, the inactive G-protein and the inactive adenylate cyclase (AC). This information is used by an algorithm that automatically draws cells on the screen: Depending on the number of signal transduction units wanted per cell, the algorithm calculates the number of cell membrane pixels needed in between the units and draws the cells accordingly (the cell membrane molecules will be extended to the left and the right of each unit in order to generate a cell with a continous border). The cell membrane molecules within the transduction unit are necessary to avoid holes in the membrane at the signal transduction sites. (b) Signal transduction unit including protein kinase C. This unit is similar to the one depicted in Fig. 4a, except that it contains an additional molecule, inactive protein kinase C (PKC). The fact that the inactive receptor is incorporated in the membrane instead of protruding does not aect the signaling. In order to make the inactive PKC accessible for modifying molecules (such as phorbolesters) both from outside and from within the cell, a background ``molecule'' was incorporated into the membrane. This allows for intracellular molecules to diuse right next to the PKC, which is necessary for an interaction.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
11
permutations in the rule table. In the above example, these are ``1,2,3,4'', ``1,2,4,3'', and so forth. The program will ®nd these 24 permutations in the rule table and then write the corresponding permutations of the product molecules ``E'', ``F'', ``G'' and ``H'', i.e. ``5,6,7,8'', ``5,6,8,7'', etc., into the 4 bytes reserved for the destination block, overwriting the default values. When the CA encounters the more common form of two molecules reacting with each other, e.g. a hormone-molecule attaching to its receptor, the program calculates the permutations by adding two other molecules to complete the four byte code (i.e. two by two array). These unspeci®ed molecules can be any of the available molecules. For example, assuming that hormone is encoded by ``3'' and receptor by ``5'', the left side of the equation ``hormone+receptor 4 active receptor+background'' will be represented by ``3, 5, 1, 1'', ``3, 5, 1, 2'', ``3, 5, 1, 3'', . . . ``3, 5, 19, 20'' and ``3, 5, 20, 20'' and sought in the rule table (the total number of combinations in this case is 4 3 202=4800). The corresponding destination blocks will then be overwritten with 4 bytes coding for one activated receptor state (coded as 6), one background state (coded as 0), and two other states corresponding to the ones calculated for the source block, e.g. ``6,0,1,1'', ``6,0,1,2'' . . . The reason for the right side of the equation to contain a distinct ``background'' is that this is a de®ned state which should not be replaced with any other molecule by the rule generating algorithm. After these four requirements are met, the user can start the simulation. An algorithm now scans the simulation window and picks up every two-by-two pixel source block in that window, arranges those 4 pixels in a row of 4 bytes, and then looks for the same 4 bytes (=source block) in the rule table. The corresponding 4 bytes of the destination block are subsequently used to replace the two by two pixels in the simulation window. The simulation algorithm also determines a random number between 0 and 32000 for every source block it encounters. Only if the reaction constant previously entered by the user is greater or equal to this number, the replacement of the source block by the destination block will take place. Scanning and replacing each block in the whole simulation window in that fashion is thought to take place simultaneously (even though it is a sequential process on non-parallel computers), and the time span of one update to the next represents the smallest unit of time in a CA, one ``generation''. All molecules, regardless of mass and volume, are represented by states of pixels in the simulation window. Moving molecules like cAMP and hormones travel in a diagonal fashion, with a change of direction after collisions with other molecules. Molecules that are thought to be not moving as much in relation to moving ones, such as membrane molecules, are ®xed in a certain position that is determined initially by the user de®ned starting pattern. The model was implemented on an Apple Macintosh Computer running system 7.x. Hyper Card (Claris, Cupertino, CA, USA) and Think C version 5.0 (Symantec, Cupertino, CA, USA) were used as programming languages.
4. Results Initial simulations demonstrated that the model comprises essential aspects of bimolecular reactions, such as hyperbolic concentration kinetics, proportionality between reaction constants and reaction kinetics, the observation of the law of mass eect and ®rst and second order
12
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
kinetics (data not shown). Another important characteristic of biologic systems, the presence of concentration gradients and their diminution over time, could be modeled by placing a multitude of molecules in a particular area of the simulation window. During the subsequent simulation, the molecules moved around the simulation window, and since there were fewer molecules to re¯ect from in the direction of the concentration gradient, the molecules on average traveled in that direction, thereby leveling the concentration gradient. If desired, membrane molecules could be de®ned in a way that promote only the passage of some molecules but not others, simulating a semipermeable membrane. Impermeable membranes Ð as in the case of the following simulations Ð are simulated by membranes that re¯ect all moving molecules. In the simulation represented by Fig. 5, a transmembranous signal transduction process is simulated. It reproduces a typical cAMP kinetic of mouse Leydig cells stimulated with CG, which is shown in Fig. 6 for comparison. As outlined in the Methods Section, the following four requirements were met before the simulation was started: The user entered the molecules, shown in Table 2, the reactions, shown in Table 3, and the starting pattern, which showed six cells with a total of 304 inactive signal transduction complexes automatically incorporated into
Fig. 5. Simulation of intracellular cAMP kinetics. Two cells with a total of 304 inactive signal-transduction units as shown in Fig. 4a are stimulated with 100 molecules of CG. Intracellular cAMP degradation is carried out by 100 phosphodiesterase (PDE) molecules, which are in turn partly inhibited by 300 molecules of PDE inhibitor. The underlying reactions as entered by the user are shown in Table 3, the rule table automatically generated based on these reactions contains 160000 individual rules (data not shown).
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
13
Fig. 6. Gonadotrophin induced intracellular cAMP-kinetics in `real' mouse-Leydig cells. Mouse Leydig cells (50000 per 0.5 ml) were incubated with two dierent concentrations of hCG (R 5 ng/ml, * 20 ng/ml) in the presence of a phosphodiestarese inhibitor (0.25 mM IBMX) for dierent periods of time at 378C. The cells were then separated from the medium by short centrifugation (1 min at 10000 g), and the determination of intracellular cAMP was carried out with a speci®c radioimmunoassay. Results are from three independent determinations; vertical bars represent 2SD, which is omitted where it is smaller than the symbol.
the cell membranes, identical to the one shown in Fig. 4a. There were also 100 CG molecules placed into the simulation window outside the cell boundaries, and 100 phosphodiesterase molecules and 300 phosphodiesterase inhibitor molecules placed within the cell boundaries. The CA then calculated the rule table with its 160000 individual rules. The simulation was then done in a window of 240 320 pixels for 750 generations, which constitutes a compromise between speed and accuracy of simulation and is therefore used in most of the simulations. The essential characteristics of the experimentally derived kinetics is captured by the simulated kinetics, such as an initially high rise in cAMP concentration, a peak value of cAMP (near to generation 350), and a slow degradation of the cAMP. This slow degradation was accomplished by introducing a phosphodiesterase and its inbitor into the model. The phosphodiesterase degrades cAMP, but it is slowed by the inhibitor. The graph shown in Fig. 5 was obtained by automatically counting the pixels with the state that represents cAMP molecules every 10 generations and then by plotting these counts against the generations. As the generations are a form of arti®cial time, such a plot resembles a time course of the concentration of whatever molecule was counted.
14
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
Table 3 Reactions that constitute the signal transduction rule used to generate the results of Fig. 5. On the left side the educts are listed, on the right side the products; reaction constants are 32000 except for reactions 2 (c = 4000) and 11 (c = 1000) Educt side of the reactions 1. hCGa+inactive receptor 2. background+active receptor 3. active receptor+inactive G-protein 4. active G-protein+inactive ACc 5. inactive receptor+active G-protein 6. active G-protein+active AC 7. active AC+cellmembrane+background+background 8. cAMPb+phosphodiesterase 9. cAMP-phosphodiesterase-complex 10. phosphodiesterase+inhibitor 11. inhibited phosphodiesterase
Product side of the reactions 4 4 4 4 4 4 4 4 4 4 4
background+active receptor hCG+inactive receptor active receptor+active G-protein active G-protein+active AC inactive receptor+inactive G-protein inactive G-protein+active AC active AC+cellmembrane+cAMP+cAMP cAMP-phosphodiesterasecomplex+background phosphodiesterase inhibited phosphodiesterase phosphodiesterase+inhibitor
a
CG: human choriogonadotrophin. cAMP: cyclic adenosine monophosphate. c AC: adenylate cyclase.
b
Fig. 7 demonstrates how the ®rst steps of signal transduction are actually simulated by the model. It is important to recognize that only a very small part of the whole simulation window is shown on an enlarged scale. Initially, all transducing molecules within the membrane are inactive. When a hormone diuses into the vicinity of a receptor, this receptor becomes activated. The hormone molecule disappears, because it is thought to be part of the active receptor (reaction no. 1 in Table 3). The active receptor now activates the G-protein, which ®nally activates the adenylate cyclase. Depending on the reaction constants, the adenylate cyclase molecules become inactive again after one or more generations. When the receptor becomes inactive, the hormone appears again and diuses away (generation x + 2 in Fig. 7). In order to incorporate ampli®cation, which is one of the principals of signal transduction, one or more of the transducing molecules will revert more slowly to the inactive state compared to how fast it became activated. Also, each active adenylate cyclase Ð at least in the rules underlying the simulation presented in this paper Ð produces two cAMP molecules per generation (generation x + 4 in Fig. 7, reaction no. 7 in Table 3). The ®nal simulations presented here focus on homologous and heterologous desensitization. This desensitization is achieved by simulating the pre-incubation of six cells either with the hormone CG (homologous, Fig. 8) or with the phorbol ester PMA (heterologous, Fig. 9). According to the reactions entered (similar to those shown in Table 3), CG activates the receptor, which activates the G-Protein, which activates adenylate cyclase as well as protein kinase C. In contrast, PMA directly activates protein kinase C. Activated protein kinase C in turn modi®es the G-protein so that it can not be activated by an active receptor (`desensitized G-protein'). The pre-incubation with CG spans about 1200 generations, the one with PMA about 600 generations. In both cases about 180 out of 240 G-proteins are desensitized after the pre-incubations, the dierent time spans being a result of dierent reactions constants involved
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
15
Fig. 7. (Caption overleaf).
16
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
in the activation of protein kinase C. Prior to re-stimulation, the simulation windows are cleared of all protein kinase C, CG and PMA molecules (simply by deleting the corresponding pixels). Re-stimulations are undertaken with either CG or AlFÿ 4 . The results depicted in Fig. 8 show clearly that pre-incubation with CG leads to a diminished response in terms of cAMPaccumulation during a second stimulation with either CG or AlFÿ 4 . A similar result is observed after pre-incubation with PMA (Fig. 9). The standard deviations were obtained by recording cAMP kinetics in three dierent simulations that contained the same cells and 100 molecules of CG. For each simulation, the CG molecules were distributed in a dierent random manner on the simulation window. 5. Discussion In this study we have described a computer model for cellular signal transduction pathways based on a cellular automaton. Although CA have been used before to model natural phenomena, they were mainly applied in the ®eld of physics, i.e. ¯uid dynamics. In recent years there were some attempts to utilize the advantages of CA in the medical ®eld, focusing on immunology, cardiology and experimental oncology. However, models to simulate signal transduction cascades have had a mathematically based approach so far, using dierential equations and special algorithms to solve them [15±17]. To our knowledge, this is the ®rst report that characterizes a CA speci®cally designed for the purpose of cellular signal transduction. The ®rst step in building a model is to clearly formulate a theory of how the dierent parts of the system to be studied interact. The CA model presented here requires the input of the molecules which play a part in the signal transduction, their concentrations, and the reaction constants of all reactions involved. Not all of this information is presently known, so some assumptions had to be made. For example, the stoichiometry of receptors to G-proteins and to adenylate cyclase molecules was arbitrary chosen to be 1:1:1 in the signal transduction units. Arbitrary numbers were also taken for other molecules, such as hormone molecules, since it was the purpose to show qualitative behavior of the cellular automaton rather than exact quantitative results. Arbitrary in this context means the absolute number of molecules of each kind taking part in the simulation, not the number of molecules that interact in a reaction, which is ®xed in the rules as exempli®ed in Table 3. The reaction constants that were used were based on our own experimental observations and
Fig. 7. Simulation of signal transduction. This sequel depicts a small part of the simulation window in the course of 5 generations during a simulation of signal transduction. A small part of the cellmembrane that contains a single signal transduction unit as explained in Fig. 5a is shown. Once a CG pixel is in direct contact with the receptor, the sourceblock that contains the receptor and the CG molecule is substituted with a destination block that contains only an activated receptor, leaving a background pixel where the CG had been. In the following generation (x + 2), this leads to activation of the G-protein, and at the same time to inactivation of the receptor, which releases the CG. In generation x + 3, the active G-proteins turns inactive again, and the adenylate cyclase becomes active. This leads to the production of two pixels of cyclic Adenosine Monophosphate (cAMP), and the inactivation of the adenylate cyclase in generation x + 4.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
17
Fig. 8. Simulation of adenylate cyclase activity after CG induced desensitization. Six cells with a total of 240 signaltransduction units as displayed in Fig. 2b were incubated for 1200 generations with (.) or without (w) 100 molecules of CG in a 184 320 pixels large simulation window. Afterwards all CG and PKC molecules were removed, and a restimulation was carried out with 100 molecules of CG (a) or AlFÿ 4 (b). The number of cAMP molecules measured are means of three dierent simulations 2SD.
published data of signaling cascades. Since the possible numbers 0 to 32000 cannot fully span all reaction constants possible in such systems, their values were assigned on a rational basis with regard to known constants and then modi®ed by trial and error. However, only a minimal amount of trial and error was necessary, which involved only the correction of reaction constants in order to have the model generate results in a reasonable time. The fact that there were no corrections to be made on the algorithms of the CA nor the fundamental reactions of
Fig. 9. Simulation of adenylate cyclase activity after PMA induced desensitization. The simulation setup was as described for Fig. 5, except that cells were ``preincubated'' with 100 molecules of phorbolester prior to removal of this PKC-stimulating agent. Restimulation was carried with CG (a) or AlFÿ 4 (b).
18
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
the model suggest that this is due to the naturalistic way CA emulate nature, i.e. generation of complex behavior based on simple rules and their application to a large number of interacting molecules. Our extension of the CA model by incorporating random numbers that in¯uence the basic reactions takes deterministic properties out of these low level calculations, whereby it most likely enhances the generation of complex behavior. A number of simpli®cations of the system were necessary to formulate the model. This was partly due to limitations of a simulation as such, and partly due to the limitations forced upon the CA by the computer hardware. Simpli®cations made in our CA include the principal restriction of all structures like cells to two dimensions, the restricted possibility of molecule movement in only four dierent directions in these two dimensions, and the restriction of only four dierent molecules interacting at one time in a given neighborhood (because it is restricted to a two by two array). In addition, characteristics such as molecular weight, molecular speed and other physical properties cannot be simulated. Furthermore, the total number of molecules being simulated is limited to 475000 at the moment, which leads to problems when very small numbers of molecules are supposed to react with each other, for example a reaction of ATP with an activated adenylate cyclase. If the simulation were to be carried out in the way that an activated adenylate cyclase would only form a cAMP molecule if it had contact with an ATP molecule, the whole simulation window would have to be ®lled with ATP to guarantee a reasonable production of cAMP. Obviously this is not possible, so simpli®cations had to be introduced to avoid these kinds of problems. In addition, due to its limitations in the number of possible reaction constants and the arbitrary selection of the ®rst reaction constant in a given set of reactions, a simple relation of the backward reaction constants and the forward reaction constants with the equilibrium constant does not exist. Nevertheless, it was demonstrated that other basic laws of chemical reactions, such as exponential reaction kinetics, the law of mass eect and the law of mass conservation are obeyed by the model. Despite these limitations, simulations were fairly accurate in the behavior that was modeled, even when only a small number of hormone molecules were used (100 molecules of CG in the simulation depicted in Figs. 8 and 9). Repeated simulations with a dierent spatial distribution of those molecules resulted in similar cAMP kinetics, as evidenced by small standard deviations of cAMP concentration at various time points. One might expect with only a few CG molecules stochastic errors of small samples, resulting in dierent cAMP kinetics. Apparently, the numbers of molecules used here were already large enough to minimize these errors. Based on experimental observations, our model ful®lls a number of important general characteristics of cellular signal transduction: . . . . . .
the hormone binds speci®cally to its receptor the hormone±receptor-complex activates a G-protein the active G-protein activates other eector molecules, such as adenylate cyclase the eector molecule produces a second messenger (in this case cAMP) there is an ampli®cation of the signal, expressed as a low hormone/cAMP ratio there is an activation of both second messenger generation and desensitization
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
19
The simulation results of the CA presented here show that CA are capable of dealing with simulations of biochemical reaction cascades, especially those involved in cellular signal transduction, in a qualitative and semi-quantitative way (semi-quantitative in the sense that reaction constants can be used to generate dierent reaction-kinetics in a network of reactions). In order to approach simulations on a true molecular scale that will predict exact concentrations of molecules at dierent time points, CA have to be further developed, mainly with respect to size and speed. Parallel to this development, knowledge about the components involved in the cellular signaling cascades will increase, enabling us to incorporate more detail into the models. Since the number of reactions and components is limited, ultimately all of these will be known. At that point the bene®t of simulations will become more obvious, because the pure knowledge of the reaction mechanisms will not be sucient to predict and thoroughly understand cellular behavior. The data presented in this paper show that CA based on the methods which we developed can serve as a basis for the development of systems that will be able to perform such simulations.
6. Summary On the basis of cellular automata (CA) models, a software speci®cally tailored to model biochemical reactions involved in cellular signal transduction was implemented on a personal computer. Although in recent years there were some attempts to utilize the advantages of CA in the biomedical ®eld, models that simulated signal transduction cascades had only utilized dierential equations and special algorithms to solve them. To our knowledge, this is the ®rst report that characterizes a cellular automaton speci®cally designed for the purpose of cellular signal transduction. Initial simulations with this model demonstrated that it comprises essential aspects of bimolecular reactions, such as hyperbolic concentration kinetics, proportionality between reaction constants and reaction kinetics, the observation of the law of mass eect and ®rst and second order kinetics. We then tailored the model for the simulation of transmembrane signal transduction in mouse Leydig cells, a process which is known to involve the hormones gondadotrophin luteotropin and/or chorionic gonadotropin. Following their interaction with cell-surface receptors, these hormones are capable of activating a transmembrane signaling process that leads to an increase of the intracellular concentration of the second messenger cyclic adenosine monophosphate (cAMP) and a subsequent stimulation of steroid-hormone biosynthesis. However, pretreatment of Leydig-cells with dierent stimulators results in a modi®cation of the signal transduction cascade leading to a diminished response of the cells during subsequent stimulations, a process known as desensitization. Recent data regarding these desensitization processes in Leydig cells were subsequently used to simulate the underlying reactions of signal transduction. The model is capable of simulating this complex behavior of the intracellular second messenger production, which involves the activation of second messenger generation and desensitization processes. Moreover, it ful®lls other important biological characteristics of cellular signal transduction, namely speci®city of hormone±receptor interactions, the interaction of signaling intermediates, and the ampli®cation of the extracellular signal, expressed as a hormone/intracellular-messenger ratio. The simulation results of the cellular automaton presented here show that CA are capable of
20
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
simulating biochemical reaction cascades, especially those involved in cellular signal transduction, in a qualitative and semi-quantitative way and suggest that quantitative simulations on a molecular level will be possible once appropriate computer hardware is available. As simulations and results of the CA are easily described and comprehended, CA are useful tools that will facilitate research in cellular signal transduction and other ®elds covering complex reaction networks. Acknowledgements The authors would like to thank Prof. Leidenberger for moral and ®nancial support to carry out this work, and Jan WuÈrthner for many programming tips. References [1] T. Tooli, N. Margolus, Cellular Automata Machines Ð A New Environment for Modeling, MIT Press, Massachusetts, 1990. [2] H. Gutowitz, Cellular Automata Ð Theory and Experiment, MIT Press, Cambridge, Massachusetts, 1991. [3] N. Margolus, T. Tooli, Cellular automata machines, in: E. Doolen (Ed.), Lattice Gas Methods for Partial Dierential Equations, SFI SISOC, Addison-Wesley, Cambridge, Massachusetts, 1990, pp. 219±249. [4] H.B.A. Sieburg, Logical dynamic systems approach to the regulation of antigen-driven lymphocyte stimulation, in: H.B.A. Sieburg (Ed.), Theoretical Immunology, Addison-Wesley, 1988, pp. 273±291. [5] H.B. Sieburg, J.A. McCutchan, O.K. Clay, L. Cabalerro, J.J. Ostlund, Simulations of HIV infection in arti®cial immune systems, Physica D 45 (1990) 208±227. [6] J.A. Reggia, S.L. Armentrout, H. Chou, Y. Peng, Simple systems that exhibit self-directed replication, Science 259 (1993) 1282±1287. [7] A.S. Qi, X. Zheng, C.Y. Du, B.S.A. An, Cellular automaton model of cancerous growth, J. Theor. Biol. 161 (1993) 1±12. [8] D.W. Frazier, P.D. Wolf, J.M. Wharton, A.S. Tang, W.M. Smith, R.E. Ideker, Stimulus induced critical point. Mechanism for electrical initiation of reentry in normal canine myocardium, Journal of Clinical Investigation 83 (1989) 1039±1052. [9] J.M. Saez, Leydig cells: endocrine, paracrine, and autocrine regulation, Endocr. Rev. 15 (1994) 574±626. [10] M. Schumacher, M. Schwarz, W. BraÈndle, Desensitation of the cAMP system in mouse Leydig cells by hCG, cholera toxin, dibutyryl cAMP: Localization of the `lesion' to the guanine nucleotide regulatory protein± adenylate cyclase complex, Mol. Cell. Endocrinol. 34 (1984) 67±80. [11] A.K. Mukhopadhyay, M. Schumacher, Inhibition of hCG-stimulated adenylate cyclase in puri®ed mouse Leydig cells by the phorbol ester PMA, FEBS Lett. 187 (1985) 56±60. [12] T. Tooli, N. Margolus, Programmable matter: Concepts and realization, Physica D 47 (1991) 263±272. [13] J.U. WuÈrthner, M. Kistler, M. Kratzmeier, A.K. Mukhopadhyay, LH/hCG-receptor is coupled to both adenylate cyclase and protein kinase C signaling pathways in isolated mouse Leydig cells, Endocrine J. 3 (1995) 579±584. [14] M. Schumacher, G. SchaÈfer, V. Lichtenberg, H. Hilz, Maximal steroidogeneic capacity of mouse leydig cells, FEBS Lett. 107 (1979) 398±402. [15] J.E.A. McIntosh, R.P. McIntosh, Mathematical Modelling and Computers in Endocrinology, Springer-Verlag, New York, 1980. [16] M.D. Houslay, The use of selective inhibtitors and computer modelling to evaluate the role of speci®c high anity cyclic AMP phosphodiesterases in the hormonal regulation of hepatocyte intracellular cyclic AMP concentratrions, Cellular Signaling 2 (1989) 85±98. [17] S. Swillens, D. Mercan, Computer simulation of a cytosolic calcium oscillator, Biochem. J. 271 (1990) 835±838.
J.U. Wurthner et al. / Computers in Biology and Medicine 30 (2000) 1±21
21
Jens U. Wurthner obtained his MD from the University of Hamburg Medical School in 1995. In the same year he began his post graduate training in internal medicine at the University of Freiburg Medical Center in the department for hematology/oncology. In 1997 he received a DFG scholarship for cancer research and subsequently joined the Laboratory of Cell Regulation and Carcinogenesis at the National Cancer Institute in Bethesda, where he is currently working on intermediates and modi®ers of the TGF-b signaling cascade. Amal Mukhopadhyay obtained his Ph.D. in Biochemistry from the Banaras H. University, Varanasi, India in 1971 and joined this University as a Lecturer of Physiology in the same year. In 1982 he moved to Hamburg to join the Institute for Hormone and Fertility Research at the University of Hamburg, Hamburg, Germany where is currently working as the Director, Division of Reproductive Endocrinology. His area of specialization is molecular basis of hormone-induced signalling mechanism in reproductive tissues. Claus J. Peimann received his diploma in physics from the University of Hamburg, Germany in 1971 and his Ph.D. in Mathematics in 1974, also University of Hamburg. He obtained his Habilitation in Medical Informatics in 1986 at the University of Hamburg Medical School and is currently an assistant professor at the Institute for Computerscience in Medicine, University of Hamburg Medical School. His special interests include simulation methods like Petri Nets and cellular automata, Hypertext/ Hypermedia and teachware.
