Quantification of Short Term Signaling by the Epidermal Growth Factor Receptor*

During the past decade, our knowledge of molecular mechanisms involved in growth factor signaling has proliferated almost explosively. However, the kinetics and control of information transfer through signaling networks remain poorly understood. This paper combines experimental kinetic analysis and computational modeling of the short term pattern of cellular responses to epidermal growth factor (EGF) in isolated hepatocytes. The experimental data show transient tyrosine phosphorylation of the EGF receptor (EGFR) and transient or sustained response patterns in multiple signaling proteins targeted by EGFR. Transient responses exhibit pronounced maxima, reached within 15–30 s of EGF stimulation and followed by a decline to relatively low (quasi-steady-state) levels. In contrast to earlier suggestions, we demonstrate that the experimentally observed transients can be accounted for without requiring receptor-mediated activation of specific tyrosine phosphatases, following EGF stimulation. The kinetic model predicts how the cellular response is controlled by the relative levels and activity states of signaling proteins and under what conditions activation patterns are transient or sustained. EGFR signaling patterns appear to be robust with respect to variations in many elemental rate constants within the range of experimentally measured values. On the other hand, we specify which changes in the kinetic scheme, rate constants, and total amounts of molecular factors involved are incompatible with the experimentally observed kinetics of signal transfer. Quantitation of signaling network responses to growth factors allows us to assess how cells process information controlling their growth and differentiation.

The epidermal growth factor receptor (EGFR) 1 belongs to the family of protein-tyrosine kinase receptors, which regulate cell growth, survival, proliferation, and differentiation (1)(2)(3). EGFR is activated by binding of epidermal growth factor (EGF) or another EGF family factor (e.g. transforming growth factor-␣).
This binding causes EGFR dimerization and rapid activation of its intrinsic tyrosine kinase followed by autophosphorylation of multiple tyrosine residues in the cytoplasmic receptor domain. Tyrosine phosphorylation of EGFR generates a biochemical message for a battery of cytoplasmic target proteins that contain characteristic protein domains, such as Src homology 2 (SH2) domains and phosphotyrosine binding domains (e.g. see . Binding and phosphorylation/activation of these proteins, e.g. growth factor receptor-binding protein 2 (Grb2), Src homology and collagen domain protein (Shc), phospholipase C-␥ (PLC␥), and others lead to a further propagation of the signal through multiple interacting pathways.
Several signaling pathways emanating from EGFR involve activation of SOS (Son of Sevenless homolog protein), the downstream target of which is Ras protein. Mitogenic signaling by EGFR is associated with Ras-dependent stimulation of mitogen-activated protein kinase cascades, leading to phosphorylation of both cytoplasmic and nuclear targets. Although a predominant role of EGFR and other tyrosine kinase receptors is stimulation of cell growth and proliferation, recent data suggest that the physiological outcome of tyrosine kinase signaling strongly depends on the timing, duration, and amplitude of activation of signaling components (2,7,8).
Initially, signaling pathways were viewed as linear relay routes, which simply transmitted and amplified signals. Now it is increasingly appreciated that signaling responses are shaped by multiple interactions of many components of signaling networks (9). A subtle difference in input signals and/or interaction kinetics may result in differential response patterns and, eventually, in alterations in gene expression by signal-regulated transcription factors. For instance, variable strength of Raf-1 activation (the first kinase of the mitogen-activated protein kinase cascade, a direct downstream target of Ras) has been linked to such opposing responses as the induction of DNA synthesis and growth inhibition (10 -13). Experiments with PC12 cells have shown that the specificity of cellular responses depends on the duration of activation of extracellular signalregulated kinase (Erk) (the terminal kinase of the mitogenactivated protein kinase cascade), e.g. whether Erk activation is transient or sustained (2, 14 -16). Therefore, signaling through the same pathway in the same cell type may result in completely different outputs depending on the amplitude and persistence of activation of signaling intermediates, i.e. on their kinetic behavior.
The kinetics (i.e. the transient and steady-state behavior) of the cellular response to EGF depends on many factors, including the number of receptors displayed on the cell surface; the concentration of the growth factor, docking, and target proteins; and their initial activity states. Moreover, other signaling pathways that share or interact with one or more components of the EGFR pathway can influence the kinetic pattern of EGFR signaling. Although a large body of data describes EGFR signaling at the molecular level, the manner in which the complex pattern of cellular responses to EGF is controlled remains poorly understood. An important reason is the lack of a quantitative description of EGFR signaling network, which hampers a careful examination of the influence of multiple factors. Detailed understanding of the dynamics of complex cellular responses requires a combination of experimental and computational approaches (17)(18)(19).
The early events of EGFR signaling, such as EGF binding and receptor autophosphorylation, binding and activation of Grb2, phosphorylation of Shc and PLC␥, and activation of SOS, develop in a time frame of seconds. There are slower processes involving receptor internalization and its subsequent degradation in lysosomes, which have an important role in EGF-induced signaling (20 -22). Activation and binding of ligands causes the recruitment of EGFR to clathrin-coated pits and transfer to endosomes. These processes are developing over time frames of minutes to hours, i.e. much more slowly than early EGFR signaling events, which evolve to a quasi-steadystate level in a time scale of seconds. Here we will study the short term response (up to 120 s) to EGF stimulation.
The aim of the present study is to give a quantitative description of the short term EGFR signaling based on a detailed kinetic scheme of the interactions of proteins and other signaling molecules involved. To this end, we combine a computational approach with experimental analysis of the time course of activation/phosphorylation of different components of the EGFR signaling pathway in isolated rat hepatocytes. We test the model against the experimental data to gain a better understanding of the factors governing the kinetics of phosphorylated signaling intermediates. In particular, the model explains why the total phosphorylated EGFR and its complexes with target proteins exhibit pronounced maxima and then descend to sustained levels, whereas the total concentration of phosphorylated forms of Shc and the concentrations of Shc-Grb2 and Shc-Grb2-SOS complexes increase monotonically, reaching a quasi-steady-state level. We demonstrate which enzyme activities and kinetic constants exert significant control over the EGFR signaling and how the transient behavior is regulated. This analysis will enable us to assess how the EGFR signaling system can process information and generate distinct outputs in response to stimuli.

Materials
Antibodies against EGFR (sheep polyclonal) and PLC␥ (mixed mouse monoclonals) were obtained from Upstate Biotechnology, Inc. (Lake Placid, NY); anti-Shc (rabbit polyclonal or mouse monoclonal), anti-Grb2 (mouse monoclonal), and anti-phosphotyrosine-horseradish peroxidase (type RC20H) were from Transduction Laboratories (Lexington, KY). Anti-phosphotyrosine-agarose conjugates (mouse monoclonal) were obtained from Sigma, and anti-IgG-horseradish peroxidase conjugates were from Pierce. Gradient gels and nitrocellulose membranes were from Bio-Rad, and detection of the Western blots was done by chemiluminescence using Supersignal reagent (Pierce). Collagenase type I was from Worthington, and bovine serum albumin fraction V and the Complete protease inhibitor mixture were obtained from Roche Molecular Biochemicals. EGF (receptor grade), protein G-Sepharose, and protein A-Sepharose were from Sigma. Other chemicals and biochemicals were obtained from Sigma or Fisher.

Cell Preparation and Incubation Conditions
Isolated hepatocytes were prepared from the livers of male Harlan Sprague Dawley rats by collagenase perfusion as described previously (23). Cell preparations were suspended in a modified Krebs-Ringer bicarbonate buffer (pH 7.4) containing NaCl (127 mM), NaHCO 3 (25 mM), KCl (4 mM), MgCl 2 (1.2 mM), potassium phosphate (1.2 mM), Hepes (10 mM, pH 7.4), CaCl 2 (1.0 mM), and glucose (15 mM) and stored on ice until use. Incubations were carried out in a shaking water bath at 37°C in capped plastic flasks in a gas phase of 95% O 2 , 5% CO 2 . Cells were preincubated at a cell density of 10 7 cells/ml in Krebs-Ringer bicarbonate buffer for 45 min to optimize receptor presentation on the cell surface, prior to stimulation with different concentrations of EGF (24). Reactions were stopped after 0, 15, 30, 45, 60, and 120 s by a 1:1 dilution of a sample of the incubation mixture with ice-cold lysis buffer containing (final concentrations) Hepes (50 mM, pH 7.5), NaCl (150 mM), EGTA (5 mM), glycerol (10%), Triton X-100 (1%), NaF (100 mM), sodium ovanadate (0.2 mM), sodium pyrophosphate (10 mM), and the commercially available protease inhibitor mixture Complete (Roche Molecular Biochemicals). After 10 min on ice, lysates were centrifuged in the cold room (4°C) in an Eppendorf microcentrifuge 5 min at top speed to remove the Triton-insoluble fraction and either used immediately or stored at Ϫ70°C until use. In some experiments, cells were lysed with a digitonin-containing buffer instead of with Triton X-100. The concentration of digitonin was 150 g/ml, equivalent to 5-7 g/mg of protein, sufficient to achieve maximal release of soluble proteins, as determined from the release of lactate dehydrogenase. Other treatment conditions were similar to those described above. The particulate fraction from digitonin-treated cells was further extracted by resuspending in lysis buffer containing 0.1% SDS plus 1% deoxycholate.

Immunoprecipitation Conditions, Gel Electrophoresis, and Western Blotting
Immunoprecipitations were carried out by a modification of the procedures described in Ref. 25. Antibody titration experiments established that maximally effective (Ͼ90%) immunoprecipitation of both unphosphorylated and phosphorylated forms of the relevant proteins was achieved at high antibody:antigen ratio. Based on these titration studies, equal volumes (25-100 l) of lysate and undiluted commercial antibody solution were mixed and incubated for a minimum of 4 h at 4°C with continuous mixing. Immune complexes with anti-EGFR antibody (sheep, polyclonal IgG) and PLC-␥ (mixed mouse monoclonal) were captured by the addition of 15 l of protein G-Sepharose added during the final hour, and anti-SHC (rabbit polyclonal IgG) complexes were captured with 15 l of protein A-Sepharose. Immune complexes were washed three times with HNTG buffer (Hepes (20 mM, pH 7.5), NaCl (150 mM), Triton X-100 (0.1%), glycerol (10%), sodium o-vanadate (0.2 mM), NaF (10 mM), and the Complete protease inhibitor mixture). Immunoprecipitates and full lysate samples were dissolved in Laemmli buffer (20% glycerol, 3% SDS, 3% ␤-mercaptoethanol, 10 mM EDTA, 0.05% bromphenol blue) and placed in a boiling water bath for 5 min. Proteins were separated by SDS-polyacrylamide gel electrophoresis on 4 -20% gradient gels electroblotted onto nitrocellulose membranes. Membranes were blocked with 1.5% bovine serum albumin in TBST (Tris 10 mM, pH 8.0, 150 mM NaCl, 0.05% Triton X-100) at room temperature. For Western blotting, all antibodies were diluted in TBST according to the manufacturer's recommendations. The membranes were incubated with primary antibody for 1 h at room temperature and then with horseradish peroxidase-conjugated secondary antibody for 30 min and washed four times for 10 min with TBST before detection by chemiluminescence.
Quantitative analysis of tyrosine phosphorylation of signaling proteins following EGF stimulation was carried out by the following procedure. Anti-phosphotyrosine (anti-Tyr(P)) immunoprecipitates of each sample were loaded onto the gels side-by-side with the corresponding samples of EGFR or PLC-␥ immunoprecipitates, after appropriate dilution with HNTG buffer to achieve a signal intensity in the same range as that of the corresponding tyrosine-phosphorylated protein on the Western blot. Alternatively, anti-Tyr(P) immunoprecipitates were loaded side by side with a sample of the corresponding total lysate (after appropriate dilution with HNTG), and the resulting nitrocellulose membranes were probed with antibodies to EGFR or SHC. This procedure allowed for quantitation of tyrosine phosphorylation of specific target proteins by normalization to the total target protein in the lysate. A similar strategy was followed to assess the extent of Grb2 coprecipitation with either EGFR or SHC proteins.
After chemiluminescence, a range of different film exposures was made for each membrane to avoid overexposure and to maintain band densities within a linear range for densitometric quantitation. Protein bands were identified according to their molecular weights and by comparison with specific immunoprecipitates. Bands were analyzed densitometrically using a Sharp JX-330 gel scanner and quantified by the Image Master ID software (Amersham Pharmacia Biotech). Results from multiple (4 -8) scans were averaged and several (three or four) immunoprecipitates from a single experiment were compared. Data are presented as the mean Ϯ S.E. for different estimates from a single experiment, representative of three or more similar experiments.

Kinetic Analysis Schematic Representation of Protein-Protein Interactions Induced by EGF Binding
For a quantitative analysis of the EGFR signaling network, an adequate description is required of the reactions that contribute to the experimentally detected protein-protein interactions and tyrosine phosphorylation events. The kinetic scheme presented in Fig. 1 forms the basis for the integration of the experimental study and the computational analysis.
In step 1, EGF binds to the extracellular domain of the monomeric EGFR (designated as R in the kinetic scheme) and forms the EGF⅐EGFR complex (designated as R a ). EGF binding drives the association of two receptor monomers into an activated receptor dimer (step 2). Recent studies (26,27) have shown that a 2:2 (EGF:EGFR) complex is the predominant form of the receptor dimer (designated as R 2 ). The phosphorylation of tyrosine residues by receptor tyrosine kinase is combined into a single step 3, yielding a form designated as RP. Although multiple tyrosine residues on the cytoplasmic tail of the receptor are targets for autophosphorylation, we did not attempt to distinguish experimentally between different phosphorylated forms of the receptor, and, as we will discuss below, the initial computational analysis also does not require a functional distinction to be made.
Step 4 is the dephosphorylation of RP, catalyzed by one or more phosphotyrosine phosphatase(s) (28,29).
Tyrosine phosphorylation triggers the binding of cytoplasmic proteins to the receptor. We consider here three proteins that directly interact with phosphotyrosine residues on the receptor, namely Grb2, Shc, and PLC␥ (4). Although several other proteins bind to the activated EGFR, it is helpful to consider a limited number of target proteins as an initial core model. It is not entirely clear whether these multiple proteins can bind simultaneously to their target phosphotyrosine residues on the same receptor molecule or whether the binding of, for example, Grb2 to the receptor hampers the binding of PLC␥ (competitive binding). The model depicted in Fig. 1 considers the binding of cytoplasmic proteins to occur by a competitive mechanism. The advantage of a model with competitive binding is that it allows us to consider receptor phosphorylation as a single step rather than monitoring different phosphorylated forms of R 2 as distinct entities. We also assume that, when Grb2, Shc, or PLC␥ are bound to EGFR, the corresponding phosphotyrosine residues are not available to receptor phosphotyrosine phosphatases. The implications of these assumptions for the dynamic pattern of EGFR signaling will be considered below. Which mechanism of interactions of EGFR and adapter proteins occurs in vivo remains to be identified.
The entire network of reactions of the receptor with its cytoplasmic target proteins can now be divided into three coupled cycles of interactions with Grb2, PLC␥, and Shc, respectively. One receptor cycle includes the binding of PLC␥ (step 5 in Fig. 1, resulting in the formation of the complex designated as R-PL) and phosphorylation of PLC␥ at two tyrosine residues by receptor tyrosine kinase (step 6, yielding the complex R-PLP). The partial cycle of the receptor is completed by the dissociation of R-PLP into phosphorylated phospholipase C␥ (PLC␥P) and RP in step 7. Tyrosine phosphorylation of PLC␥ is thought to be necessary for its activation and the subsequent formation of inositol 1,4,5-trisphosphate and generation of a Ca 2ϩ response (30,31). PLC␥P can translocate to cytoskeletal or membrane structures (step 25), which yields bound PLC␥P-I (32,33).
Another partial receptor cycle starts with the binding of Grb2 to a receptor phosphotyrosine (step 9, forming the complex R-G). The complex of the EGF receptor with the adapter protein Grb2 is a branch point that leads to several signaling pathways through binding to different potential targets. Here we consider the link of EGFR to the Ras signaling pathway. The SH3 domains of Grb2 bind to proline-rich regions of the Ras-specific GDP-GTP exchange factor SOS. In step 10, SOS binds to the receptor-bound Grb2, resulting in the formation of the ternary complex R-G-S. The binding of SOS to the EGFR-Grb2 complex localizes SOS in the vicinity of Ras, which is anchored to the cell membrane. The ternary complex R-G-S dissociates (step 11), yielding the phosphorylated receptor (RP) and the complex G-S, which further dissociates into Grb2 and SOS (step 12).
The final EGFR cycle considered here includes the formation of the complex of Shc with EGFR (R-SH) (step 13 in Fig. 1) and its subsequent phosphorylation at Tyr 317 by receptor tyrosine kinase (step 14, yielding R-ShP). This allows Grb2 to also bind to EGFR indirectly through phosphorylated Shc, forming a ternary complex (R-Sh-G) (step 17). There are three embedded EGFR cycles that involve Shc protein. The shortest of these cycles is completed in step 15, where the complex R-ShP dissociates, yielding the phosphorylated receptor (RP) and phosphorylated Shc (ShP). The second cycle is completed in step 18, where the ternary complex R-Sh-G dissociates into RP and the complex Sh-G. The longest of the three embedded cycles includes SOS binding to R-Sh-G, leading to the formation of a four-protein complex, R-Sh-G-S (step 19). The complex R-Sh-G-S can also be formed by association of R-ShP and G-S complexes in step 24. The third cycle is completed in step 20, where the complex R-Sh-G-S dissociates, releasing the phosphorylated receptor (RP) and the complex Sh-G-S.
It is unknown whether the binding of the phosphorylated target proteins to EGFR protects them against specific phosphatases. The kinetic scheme of Fig. 1 assumes that PLC␥P and ShP are dephosphorylated only after they dissociate from the receptor (steps 16 and 8). However, this assumption is not critical, provided the dephosphorylation of bound target proteins proceeds no faster than that of their unbound phosphorylated forms.
After phosphorylated Shc dissociates from the receptor (ShP), it retains its ability to bind various SH2 domain-containing targets. The remaining steps in Fig. 1 constitute the cycle of ShP. The scheme shows that Grb2 binds to ShP, forming the complex Sh-G (step 21). The GDP-GTP exchange factor SOS is able to bind to Grb2 complexed with phosphorylated Shc, forming the ternary complex Sh-G-S (step 22). The dissociation of the complex Sh-G-S yields G-S and ShP (step 23).

Derivation of a Kinetic Model
Kinetic Equations-In order to integrate the experimental observations in a description of the dynamic behavior of the EGFR signaling network, we converted the reaction scheme of Fig. 1 into a set of mathematical equations known as chemical kinetics equations (34). For changes with time of the concentration of any component, e.g. the receptor form RP, one can write the following.
Rate of change of RP concentration ϭ total rate of RP production Ϫ total rate of RP consumption (Eq. 1) Here the total rate is the sum of the rates that produce or consume RP according to the kinetic diagram. For instance, the total rate of RP production equals the sum of the (net) rates of six steps (steps 3, 7, 11, 15, 18, and 20; see  Table I. Kinetic equations are usually written in terms of concentrations (not of mole numbers), since the reaction rates are functions of concentrations. If the same compound participates in reactions taking place in different compartments with different volumes, the effective concentration of that compound will be different depending on the volume of the corresponding compartment.
Step 1 (EGF binding to EGFR) could be considered as taking place in the extracellular compartment with a given initial concentration of EGF. The concentration of EGFR in the extracellular compartment would then be calculated as the number of the receptors on the cell surface divided by the (average) volume of incubation medium per cell (V m ). In step 2, association and dissociation of the receptor monomers occurs in the cell membrane. All other steps are considered as taking place in the cytosolic compartment. Therefore, the same mole number of EGFR would give rise to three EGFR concentrations (representing the different compartments). However, for computational purposes, it is more convenient to deal only with a single concentration of EGFR related to the cytoplasmic water volume (V cw ) of the cell. This requires rescaling the rate constants of steps 1 and 2. For the purpose of this rescaling, the EGF concentration in the model was also related to the cytoplasmic water volume; i.e. [EGF] in the experimental medium was multiplied by the ratio V m /V cw (see Table II). Typically, there were 10 7 cells/ml in our experiments (see "Cell Preparation and Incubation Conditions"); therefore, V m ϭ 10 Ϫ7 ml. Assuming the diameter of a hepatocyte of 20 m and a cytoplasmic water volume of about 70% of total intracellular volume, V m /V cw ϭ 33.3.
Conserved Moieties-In the reaction network described by the equations listed in Table I, the EGF moiety and the protein moieties are conserved. This assumption is justified for the short term responses considered here. Let [EGFR] total be the total concentrations of EGFR forms. Then the following is true.
Assuming that 60 -80% out of the total of 1-3⅐10 5 EGF receptors/cell (35)(36)(37) is displayed on the cell membrane, the total concentration of surface-expressed EGFR, translated to the cytoplasm water volume, is about 100 nM. Five other moieties conserved in the EGFR signaling reactions include the total concentrations of PLC␥, Grb2, Shc, and SOS proteins and EGF, designated below by [ (Table II).
Thermodynamic Restrictions along Cyclic Pathways in the Kinetic Scheme-If a kinetic scheme includes "true" cycles, in which the initial and final states are identical, the equilibrium constants of the reactions along any cycle satisfy so-called "detailed balance" relationships (e.g. see Refs. 38 and 39). These detailed balance relations require the product of the equilibrium constants along a cycle to be equal to 1, since at equilibrium the net flux through any cycle vanishes. Therefore, such relations decrease the number of independent rate constants in a kinetic model. The kinetic scheme in Fig. 1 demonstrates that the progression along steps 9 -12 (in the positive direction) completes a cycle without any concomitant transformations and changes in the free energy. Hence, the following restriction exists on the kinetic constants.
Further examination of the kinetic scheme in Fig. 1 shows additional reaction cycles that imply the following constraints.
a In order to rescale step 1 to the cytoplasmic water volume (V cw ), k 1 was multiplied by the factor V cw /V m ϭ 0.03, equal to the ratio of the cytoplasmic water volume and the volume of incubation medium per cell (V m ).
values of 0.4 nM (37) and of 0.03 and 0.29 nM for high and low affinity sites (36), respectively, have been reported. Our recent data also demonstrate that in intact hepatocytes, EGFR autophosphorylation saturates at EGF concentrations of 5-10 nM (25), whereas in Triton-solubilized cells maximal activation of EGFR tyrosine kinase activity requires 500 -1000 nM EGF. These findings indicate that the K d for EGF in intact hepatocytes should be well below the value of 100 -500 nM measured for the solubilized receptor. In the kinetic model, we have used K d ϭ 0.6 nM for EGF binding to intact liver cells. This represents an average value of literature data on binding studies in hepatocytes and is compatible with our experimental data on EGFR autophosphorylation.
Whereas the K d values are important for the (quasi)equilibrium conditions, a knowledge of the rate constants of the forward and backward reactions is required to describe the temporal (and steady-state) behavior. The association and dissociation steps are characterized by second-order and first-order rate constants, respectively. For EGF binding to the recombinant soluble extracellular binding domain of EGF receptor, the "on" (association) and "off" (dissociation) rate constants were reported to be k 1 ϭ 1.5⅐10 Ϫ4 nM Ϫ1 s Ϫ1 and k -1 ϭ 0.06 s Ϫ1 , and the K d ϭ k -1 /k 1 ϭ 400 nM (43). The K d value of 0.6 nM, characteristic for membrane-bound receptor can be obtained if the reported (43) magnitude of k 1 is increased or k -1 is decreased by a factor of about 600. The characteristic (relaxation) time of the EGF binding reaction is 1/(k -1 ϩ k 1 ⅐[EGF]). A decrease in k -1 by a factor of 600 leads to a relaxation time equal to about 2500 s (for 2 nM EGF), which is about 3 orders of magnitude higher than the characteristic time of experimentally observed responses of the entire signaling network (see Figs. 2 and 3). Therefore, we conclude that the value of the off rate constant, k -1 , is unlikely to be much less than the value reported in Ref. 43. On the other hand, the on rate constant, k 1 , could be substantially higher for the receptor in situ because of the decreased orientation restrictions on the positions of encountering molecules. Indeed, the values of this constant determined in intact fetal rat lung cells (44) and in human fibroblasts (45) were 20 and 35 times greater than the value reported by Zhou et al. (43). Taking k -1 ϭ 0.06 s Ϫ1 (43,45) and K d ϭ 0.6 nM gives k 1 ϭ 0.1 nM Ϫ1 s Ϫ1 . With these values of the rate constants, the characteristic time of the binding reaction with 2 nM EGF is less than 4 s.
Receptor Dimerization-Aggregation of activated receptor monomers (R a ) into a dimer (R 2 ) is brought about by the random lateral diffusion of R a in the cell membrane. The lateral diffusion coefficient (D R ) reported for EGFR is about 1-2⅐10 Ϫ10 cm 2 s Ϫ1 (46), in line with the values determined for various membrane proteins, which are typically in the range of 5⅐10 Ϫ9 to 10 Ϫ10 cm 2 s Ϫ1 (47)(48)(49). Substituting the reported value of D R into equations for diffusion-limited rate constants in two dimensions (50 -53) and relating the collision rate to a unit of cytoplasmic water volume, we calculated the diffusion limit for the second-order rate constant k 2 to be 1-0.02 nM Ϫ1 s Ϫ1 . Assuming that the dimerization rate is below the lower limit of the encounter rate, we have taken k 2 ϭ 0.01 nM Ϫ1 s Ϫ1 . Given a K d of EGFR dimerization of 10 nM, k -2 ϭ 0.1 s Ϫ1 .
Rate Constants of Phosphorylation and Protein Binding-In a living cell, the ATP concentration is much higher than the Michaelis constant of the receptor kinase for ATP (54,55). Therefore, the rate of tyrosine phosphorylation of the receptor (as in step 3) or bound target proteins (as in steps 6 and 14) can be kinetically characterized by pseudo-firstorder rate constants. Importantly, the standard free energy differences of the tyrosine phosphorylation reactions are low, so that the equilibrium constants are of the order of unity (56 -58). Consequently, the phosphorylation steps catalyzed by EGF receptor kinase are considered reversible (Table II), and the effective rate constants will depend on the ATP:ADP ratio, which is assumed constant. By contrast, the phosphatase reactions (steps 4, 16,8) can be considered as (kinetically) irreversible. The phosphatases are assumed to follow Michaelis-Menten kinetics (29,59), and the concentration of inorganic phosphate is considered constant.
The association of protein molecules into dimers or larger complexes occurs with typical rate constants on the order of 10 Ϫ4 to 10 Ϫ1 nM Ϫ1 s Ϫ1 (60 -62). The rate constants reported for the association of the p85 subunit of phosphatidylinositol 3-kinase with two different phosphopeptides, corresponding to phosphotyrosine sites of the plateletderived growth factor ␤-receptor were 1.9⅐10 Ϫ3 and 9.2⅐10 Ϫ3 nM Ϫ1 s Ϫ1 (63). The binding of SH2 domains of the p85 subunit to phosphotyrosine sequences derived from the insulin receptor substrate-1 was characterized by association rate constants of 3⅐10 Ϫ2 to 4⅐10 Ϫ1 nM Ϫ1 s Ϫ1 for two different phosphopeptides and N-terminal and C-terminal SH2 domains of p85 (64). The dissociation rate constants were observed to be 0.1 s Ϫ1 for platelet-derived growth factor ␤-receptor-derived peptides (63) and 0.1-0.2 s Ϫ1 for insulin receptor substrate-1-derived sequences (64). The dissociation equilibrium constants appeared to be 14 -50 nM for platelet-derived growth factor-derived peptides (63) and 0.3-3 nM for insulin receptor substrate-1 peptides (64). In a recent study (65), a much higher K d of 300 nM for the interaction of a platelet-derived growth factor ␤-receptor-derived peptide with the N-terminal SH2 domain of the p85 subunit of phosphatidylinositol 3-kinase has been reported. The discrepancy between these and some other literature data (66 -69) regarding the binding affinities of the SH2 domains can be explained by the differences in the experimental techniques, the SH2 domains, and the phosphopeptide sequences studied (65). Since data for on and off rate constants for the EGFR interactions with its target proteins in situ are unavailable, the corresponding rate constants were assumed to be in the same range as those reported for the binding of SH2 domains to phosphopeptides (see Table II).
Grb2 interacts with SOS (step 12) with a high affinity through the N-terminal SH3 domain (70). The off rate constant of Grb2-SOS complex is orders of magnitude slower than the off rate constants for the interactions of SH2 domains with phosphotyrosine peptides (70). Fast dissociation rates at the phosphorylated receptor sites are important for rapid exchange of ligands (63,64). It has been reported that the Grb2-SOS complex binds to both EGFR-and Shc-derived phosphopeptides with higher affinity than Grb2 alone (71). These experimental data restrict the K d values of the corresponding reactions in the kinetic scheme ( Fig. 1) Fig. 1 are collected in Table II.

Experimental Analysis of EGFR Signaling
The time course of the cellular response to EGF was followed in freshly isolated hepatocytes by measuring tyrosine phosphorylation and protein-protein interactions of signaling intermediates described in Fig. 1 after stimulation with different EGF concentrations (20, 2, or 0.2 nM) for 15, 30, 45, 60, and 120 s. EGFR phosphorylation was determined by two different approaches. First, EGFR protein was immunoprecipitated using an anti-EGFR antibody that recognized both phosphorylated and unphosphorylated receptor, and EGFR phosphorylation was determined by probing Western blots with anti-Tyr(P) antibody (25). Alternatively, tyrosine-phosphorylated proteins were immunoprecipitated using an anti-Tyr(P) antibody, and membranes were probed for EGFR protein by Western blotting using an anti-EGFR antibody (or, for other phosphorylated proteins, using appropriate antibodies). In some experiments, EGFR protein in the anti-Tyr(P) immunoprecipitates was compared with EGFR protein in samples of the total lysate run in parallel on the same membranes, or the anti-Tyr(P) immunoprecipitates were compared on the same membranes with the anti-EGFR immunoprecipitates, using the same anti-EGFR antibody for detection. These experiments provided estimates of the phosphorylated EGFR protein as a fraction of total EGFR protein in the lysate. Fig. 2 shows total phosphorylated EGFR as a fraction of the total EGFR protein in the lysate at different times after stimulation with EGF. The kinetic scheme in Fig. 1 indicates that the following EGFR forms contributed to the bands of the phosphorylated receptor.
Activation of hepatocytes with a saturating concentration of EGF (20 nM) elicited a rapid response of receptor autophosphorylation. The peak EGFR phosphorylation level was reached within 15 s and was equivalent to 50 -70% of the total detectable receptor protein in the cells. It declined to reach a quasisteady-state phosphorylation level of 15-20% of the total receptor population by 2 min. This sustained level was maintained during further incubation up to 30 min (not shown). A lower EGF concentration (2 nM) also induced a tran-sient receptor phosphorylation response, but the peak level was significantly less (35% of total EGFR). A much lower peak phosphorylation level (less than 10%) was obtained at 0.2 nM EGF. The peak phosphorylation level detected in anti-EGFR immunoprecipitates and in anti-Tyr(P) immunoprecipitates after a 15-s stimulation with EGF was not significantly different, indicating that conditions used for immunoprecipitation were equally effective with either antibody and that the data describing early events in EGFR signaling were not significantly affected by the method of measurement. Fig. 3 shows the time course of the EGF (20 nM)-induced activation of several downstream signaling events, as detected by tyrosine phosphorylation of PLC␥ and Shc proteins (A and B), and by Grb2 coprecipitation with EGFR and Shc (C). Tyrosine phosphorylation of PLC␥ (Fig. 3A) was measured by immunoprecipitation with anti-Tyr(P) antibody and detection by Western blotting with anti-PLC␥ antibody. The quantity of PLC␥ protein detected in these immunoprecipitates was compared with total PLC␥ protein obtained after immunoprecipitation with anti-PLC␥ antibody and run in parallel on the same gels. The PLC␥ band in the anti-Tyr(P) immunoprecipitates includes the following complexes (see Fig. 1).
Total phosphorylated PLC␥ ϭ ͓R-PLP͔ ϩ ͓PLC␥P͔ (Eq. 14) The time course of tyrosine phosphorylation of Shc (Fig. 3B) was also measured in anti-Tyr(P) immunoprecipitates using anti-Shc antibodies for detection. The predominant Shc isoforms detected in the liver cell lysate include a 46-and a 52-kDa form. Both isoforms become tyrosine-phosphorylated in response to EGF with approximately similar kinetics (72). The density of the corresponding bands was compared with the Shc protein bands detected in the total lysate analyzed in parallel on the same gel. The phosphorylated Shc protein detected in these analyses reflects the sum of the following concentrations (see Fig. 1).
Grb2 coprecipitation with EGFR and with Shc ( Fig. 3C) was measured by immunoprecipitating cell lysates with anti-EGFR antibody and with anti-Shc antibody and detecting coprecipitated Grb2 by Western blotting with Grb2 antibody. According to the kinetic scheme of Fig. 1, the corresponding bands in the gels include the following complexes.
Total concentration of Grb2 bound to EGFR forms ϭ ͓R-G͔ ϩ ͓R-G-S͔ Quantification of the bands was done by comparison with the total cell lysates analyzed in parallel on the same gels. In agreement with our earlier findings (25), a larger fraction of Grb2 was bound to Shc than to EGFR (Fig. 3C). A striking feature of the early responses to EGF is the pronounced maximum in the concentrations of phosphorylated EGFR, Grb2 coprecipitated with EGFR, and phosphorylated PLC␥ (see Figs. 2 and 3). Interestingly, the time course of EGFR phosphorylation correlates with the time course of Grb2 bound to EGFR and of phosphorylated PLC␥, suggesting that the events occurring in the different branches of the kinetic scheme of Fig. 1 were partially synchronized. The observed pattern of early EGFR signaling raises several questions, such as the following. Why do tyrosine phosphorylated EGFR and PLC␥, as well as the total concentration of Grb2-EGFR complexes, exhibit pronounced peaks and then descend to rela- tively low sustained levels despite continuous EGF stimulation? At the same time, why does the total concentration of phosphorylated forms of Shc and of Grb2-Shc complexes increase monotonically, reaching a quasistationary level? Why does the amount of Grb2 bound to Shc significantly exceed that of Grb2 bound to EGFR? Computational kinetic analysis provides a tool to answer these questions.

Computational Kinetic Analysis of EGFR Signaling Time Course of Receptor Phosphorylation and Target Protein
Recruitment-The time course of responses to EGF was computed and compared with the experimental observations. Fig.  4A (solid lines) illustrates how the total concentrations of phosphorylated receptor forms depend on the duration of hepatocyte stimulation by 20, 5, and 2 nM EGF. These transients demonstrate a good fit to the experimentally observed time course (Fig. 2B), exhibiting a marked decline in the total phosphorylated EGFR following an initial peak. Taking into account that there were 10 7 cells in 1 ml of incubation medium, computational analysis confirmed that 20 nM EGF in the medium is a saturating concentration for EGFR signaling (cf. lines 1 and 2 in Fig. 4A). The peak level of total phosphorylated EGFR, normalized to cytoplasmic water space was about 80 nM at saturating EGF concentration (i.e. 80% of the surface-expressed EGFR, corresponding to about 50% of the total EGFR population) and 50 nM with 2 nM EGF. Previously, in order to explain the early peaks in the transients, the rapid burst of tyrosine phosphorylation of EGFR and/or some cytosolic targets was assumed to cause the activation of tyrosine phosphatases that dephosphorylate Tyr(P) residues (e.g. Ref. 25). Importantly, the kinetic model indicates that this assumption is not required if binding of a target molecule to a Tyr(P) residue of EGFR protects the residue against a constitutive phosphatase activity. During the time interval when the phosphorylated EGFR proceeds through its duty cycles (steps 5-20 in Fig. 1), Tyr(P) residues occupied by their ligand proteins are protected against dephosphorylation. The total concentration of these receptor forms begins to significantly exceed the concentration of the ligand-free form, RP (which is rapidly dephosphorylated), resulting in an effective increase in tyrosine phosphorylation of the receptor. The completion of receptor cycles returns the receptor to the ligand-free RP form, hence increasing the RP concentration relative to that of nonphosphorylated dimer R 2 . Since the phosphatase(s) continue to dephosphorylate RP, the dephosphorylation rate (Fig. 4B, line 2) begins to exceed the rate of R 2 phosphorylation by tyrosine kinase (line 1), and the level of tyrosine phosphorylation of EGFR decreases. By contrast, when we assumed that Tyr(P) residues occupied by ligands are still accessible to phosphatase (which, therefore, effectively competes with the ligands), the experimentally observed maxima did not appear in the simulated responses to EGF (Fig. 4A, dashed line).
The total concentrations of phosphorylated Shc and of Grb2 coprecipitated with phosphorylated Shc do not exhibit a marked maximum (Fig. 5A), and they reach a quasistationary level, in agreement with our experimental observations (Fig.  3B). Importantly, the kinetic model explains why these transients differ so markedly from the transients of the total phosphorylated EGFR (Fig. 4A) and Grb2 coprecipitated with EGFR (Fig. 5B). Computations show that the total phosphorylated Shc bound to EGFR (i.e. [R-ShP] ϩ [R-Sh-G] ϩ [R-Sh-G-S]) exhibits a pronounced peak, descending then to a low sustained level (Fig. 5B). We conclude that the almost monotonic increase in the total phosphorylated Shc is brought about by the accumulation of the phosphorylated forms dissociated from the receptor, i.e.
The progress of phosphorylated receptor through its cycles can be described in terms of a wave propagation. Indeed, the transients of the concentration of phosphorylated EGFR forms and of the target proteins bound to EGFR behave as a single wave. A decrease in free (nonactivated) forms of the target proteins after EGF stimulation prevents the repetition of such waves, driven by tyrosine phosphorylation of the receptor at the expense of ATP hydrolysis. Because the computational model does not include the process of EGFR internalization, the completion of this transient process leads to steady-state signaling. Computations show that the quasistationary levels reached by about 2 or 3 min are maintained during the following 30 min, and this is confirmed by experimental observations (25).
It is believed that a key component of the activation of SOS is its recruitment by EGFR to the plasma membrane, where Ras protein is located (73)(74)(75). The binding of SOS to EGFR is mediated by Grb2, which forms a stable complex with SOS (step 12) even in the absence of EGF stimulation (70). Since the reported K d for this interaction is in the nanomolar range (70), a substantial fraction of SOS (more than 50%) should exist in complex with Grb2 in resting cells (G-S in Fig. 1). The kinetic model shows that after EGF stimulation the concentration of the free complex G-S decreases, as G-S binds to phosphorylated receptor. The total concentrations of SOS bound to EGFR (R-G-S and R-Sh-G-S) and to the phosphorylated Shc (Sh-G-S) exhibit transient and monotonic increases, respectively (Fig.  5C). The interaction of SOS with Ras may be more effective (in terms of the formation of the productive complex) when SOS is brought in close vicinity to the membrane-bound Ras protein than it would have been when it depends on collisions with Ras from the cytosol (76). Hence, the transient response of SOS complexed with EGFR (Fig. 5C) may result in a transient activation of Ras. It is also possible that SOS can be targeted (through Shc) to other scaffolding proteins to generate additional Ras activation signals, which can be separately controlled.
Using the kinetic model, it is instructive to monitor how the rates of individual steps change with time. The rates of steps of the receptor cycles involving EGFR interaction with PLC␥ and Shc increase to peak values and then decrease to low (less than 1 nM/s) or close to zero sustained (stationary) values. Remarkably, the rates of phosphorylation of the receptor (by EGFR intrinsic tyrosine kinase in step 3) and its dephosphorylation by phosphotyrosine phosphatase(s) (step 4) do not decrease to zero with time (Fig. 4B). On the contrary, they reach rather high sustained values. The phosphorylation and dephosphorylation cycle involving steps 3 and 4 is an ATP consumer. Our computations showed that during a sustained EGFR signaling, the energy demand is less than 0.1% of the total ATP production in hepatocytes. The stationary rates of dephosphorylation of ShP (step 16) or PLC␥P (step 8) are relatively low (less than 1 nM/s).
Origin of the Transient Response of PLC␥-The computa-tional model indicates that the experimentally observed rapid transient phosphorylation of PLC␥ (Fig. 3A) would require some form of deactivation of phosphorylated PLC␥ after its dissociation from EGFR in step 7 or the disabling of recurrent phosphorylation of PLC␥ after its dephosphorylation in step 8 (Fig. 1). If PLC␥P concentration would accumulate and increase significantly above the R-PLP concentration, the time course of accumulation of total phosphorylated PLC␥ should be monotonic (similar to that of total phosphorylated ShP). Computations show that the dynamics of the PLC␥ cycle does not fit the experimentally observed transients, unless it is assumed that PLC␥P undergoes an additional transformation preventing the immediate conversion to free unphosphorylated PLC␥. Binding of PLC␥P to the cytoskeleton or membranes, or to any other structural component of the cell can function as such a transformation. Indeed, binding of PLC␥P to the cytoskeleton The dashed line corresponds to the assumption that binding of ligands does not protect Tyr(P) residues of EGFR against the phosphatase. B, 1, rate of EGFR autophosphorylation (step 3 in Fig. 1), 2, dephosphorylation rate (step 4).

FIG. 5. Computation of the time course of downstream EGF signaling in hepatocytes.
A, total phosphorylated Shc (lines 1 and 2) and total Grb2 coprecipitated with Shc (lines 3 and 4). B, total phosphorylated Shc bound to EGFR (lines 1 and 2) and total Grb2 bound to EGFR (lines 3 and 4). C, total (activated) SOS bound to EGFR (lines 1 and 2)  has been proposed in experimental studies carried out on hepatocytes, maintained in culture for 24 h (32,33). Another possibility is that after its dephosphorylation in step 8, PLC␥ would no longer be available for interactions with EGFR, which impedes its further tyrosine phosphorylation by the receptor kinase. Fig. 5D compares the time course of the total phosphorylated PLC␥ when we assume a translocation (step 25) of PLC␥P to a structural element of the cell (solid lines) and in the absence of such a process (dashed line).
We tested experimentally to what extent binding of PLC␥ to cellular constituents could account for this response pattern. In the experiment of Fig. 6, isolated hepatocytes were stimulated with EGF (20 nM) for periods ranging from 15 s to 60 min, and cells were permeabilized with digitonin (150 g/ml, equivalent to approximately 7 g/mg of protein). At this concentration, digitonin selectively permeabilizes the plasma membrane and causes release of soluble proteins from the cytoplasm, while bound or compartmentalized proteins are retained in the particulate fraction. Specifically, all EGFR protein is recovered from the particulate fraction, indicating that digitonin treatment caused no significant solubilization of plasma membrane proteins.
The analysis of PLC␥ in the soluble and particulate fraction of digitonin-treated hepatocytes is shown in Fig. 6. The data of Fig. 6 demonstrate that the vast excess (Ͼ90%) of PLC␥ is available free in the cytosol in the dephosphorylated form, yet is excluded from sustained EGFR-mediated phosphorylation. This finding suggests that other modifications of PLC␥ or EGFR occur following EGF stimulation of the cells that prevent their interaction. However, analysis of PLC␥ with phosphoserine or phosphothreonine antibodies detected only a modest level of serine/threonine phosphorylation that was entirely confined to the PLC␥ bound to the particulate fraction. Another possible mechanism is suggested by a recent report (77) that demonstrates that phosphatidylinositol 1,4,5-trisphosphate, the product of phosphatidylinositol 3-kinase, can bind to the PLC␥ SH2 domain and inhibit its binding to phosphotyrosine residues on growth factor receptors. This mechanism may be involved in localizing PLC␥ in the vicinity of its substrate and also result in suppressing PLC␥ availability for binding to the phosphorylated growth factor receptor. To what extent this mechanism contributes to restricting access of PLC␥ to the activated EGFR in hepatocytes is currently under investigation.

Sensitivity of the Dynamic Pattern to Variations in Rate
Constants-The dynamics of the EGFR signaling appears to be robust to significant changes in the rate constants of the protein interactions involved (cf. Ref. 78). Typically, a severalfold (in many cases 1 or even 2 orders of magnitude) variation of a rate constant does not result in significant changes of the response to EGF. However, not all of the rate constants can be arbitrarily changed, and certainly, a simultaneous alteration of all of the rate constants does result in a marked change of the response dynamics. To come to grips with the latter issue, we return to the equations in Table I that describe the time course of signal propagation. The time derivative of the concentration of a component enters the left-hand side of these equations, and each term on the right-hand side is the rate of production or consumption of that component in a particular process. A simultaneous, say 2-fold, change in all of the rate constants results in multiplication of the right-hand side of every equation in Table I by a factor of 2 (note that V max of the phosphatases also increase twice after a 2-fold increase in all of the elemental rate constants, whereas Michaelis constants do not change). This multiplication is equivalent to scaling of the time, and exactly the same levels of cellular tyrosine phosphorylations will be now reached twice as fast. An important feature of the kinetic behavior of EGFR signaling results from the bimolecular nature of protein-protein interactions. A simultaneous, e.g. n-fold increase in the amounts of the proteins involved in EGFR signaling and an n-fold decrease in the second-order rate constants (of protein association reactions) at unchanged values of all the first-order rate constants will not change the time course of the signal propagation (during this procedure, the phosphatase reactions should be also considered at the level of elemental processes (79)).
Dependence of the Responses to EGF on Relative Abundance of Signaling Proteins-The kinetic model emphasizes that the dynamic pattern of signal propagation strongly depends on the relative abundance of molecular factors involved in the EGFR pathway. Because the cellular concentrations of signaling molecules such as Shc, Grb2, and SOS have not yet been determined (70), we employed a computational analysis to determine the concentration range over which the calculated responses to EGF are consistent with our experimental observations. The time course of phosphorylation/activation responses to EGF appeared to be more sensitive to variations in relative concentrations of signaling proteins than to most variations of the kinetic constants. This can be explained by a competition of various adapter/target proteins for EGFR and by nonlinear interactions leading to the formation of multiprotein complexes. The model suggests that in hepatocytes, the total concentrations of Shc and Grb2 do not differ by more than 1 or 2 orders of magnitude from the total EGFR concentration (rescaled to cell water volume), and [Shc] total and [Grb2] total exceed [SOS] total (Table II). As yet, there are no solid experimental data to assess the validity of this conclusion. Fig. 7 illustrates how the cellular content of Shc and Grb2 affects the time pattern of EGFR-Grb2 coprecipitation. Fig. 7 (line 1) shows that a 4-fold decrease in [Shc] total completely eliminates the peak in the time course of the total Grb2 bound to EGFR, in sharp contrast with experimental observations (compare Fig. 7 with Fig. 3C). Interestingly, similar changes in the response pattern to EGF are brought about by a 4-fold increase in [Grb2] total (Fig. 7, line 2). These results suggest that the Shc:Grb2 ratio is an important controlling factor of the EGFR signaling response.
Using the kinetic model, we examined how the relative abundance of EGFR affects the time course of cellular response to EGF. An increase in EGFR content by factors ranging from 2 FIG. 6. Distribution and tyrosine phosphorylation of PLC␥ after EGF stimulation. Isolated hepatocytes were stimulated with EGF (20 nM) for periods of 15 s to 60 min and permeabilized with digitonin. Supernatants were removed (soluble fraction), and pellets were reextracted with lysis buffer containing 0.1% SDS, 1% deoxycholate (particulate fraction). Immunoprecipitation of PLC␥ or tyrosine-phosphorylated proteins was carried out with antibodies against PLC␥ (A) and phosphotyrosine (B), respectively, as described under "Experimental Procedures" and analyzed by Western blotting using anti-PLC␥ antibodies (anti-PLC␥ Ab).
up to 10 did not change the qualitative behavior of the response. However, there were significant quantitative changes. For Shc phosphorylation, overexpression of EGFR markedly shortened the period of time required to reach a high sustained level. For EGFR and PLC␥ phosphorylation and for the concentration of Grb2-EGFR complexes, an increase in the amount of active receptors decreased the time required to descend from a high peak to a relatively low sustained level. The maximal levels of phosphorylation of EGFR and PLC␥ and of total Grb2 bound to EGFR rose with an increase in the total EGFR level. The peak:sustained concentration ratio decreased for phosphorylated EGFR and for Grb2 bound to EGFR but increased for phosphorylated PLC␥. A significant decrease in EGFR availability brought about both quantitative and qualitative changes in response behavior. First, it substantially changed the time frame of the progress of EGFR signaling. For instance, the time of relaxation of EGFR autophosphorylation to a sustained level increased from about 1 min to over 20 min with a 90% decrease in active EGFR content. Second, for all transients, the ratio of peak to stationary value sharply decreased with inactivation of EGFR, and a 90% decrease in [EGFR] total eliminates the transient peak in phosphorylated PLC␥. These findings illustrate the functional importance of short term and long term variations in EGFR availability (by changes in surface expression or receptor inactivation or down-regulation) for the dynamic response patterns in the cell. DISCUSSION Protein phosphorylation that results from EGF stimulation of cellular receptors is a dynamic process. It reflects not only activation of the receptor protein kinase but also the interactions between different signaling components and the activities of various phosphatases. Downstream EGFR signaling depends on how the phosphorylation of adapter and target molecules develops over time, and both transient and sustained levels of this phosphorylation/activation ultimately depend on the entire network of signaling reactions. However, there is a gap between our knowledge of signaling events at the molecular level and our understanding of how the cellular response is integrated to achieve the desired physiological outcome. In order to fill this gap, we have combined experimental analysis with a computational kinetic approach, with the goal of creating a unifying framework for studying receptor tyrosine kinase signaling.
Experimental analysis of the time course of the response to EGF stimulation in rat hepatocytes demonstrated a rapid burst in receptor phosphorylation and accumulation of phosphorylated/activated target proteins, which occurs as early as within 15-30 s following EGF stimulation (Figs. 2 and 3). The time resolution of the experiments was not sufficient to show the initial rate of increase in EGFR phosphorylation or activation of its target proteins. A more gradual increase can be seen at temperatures lower than 37°C used in these experiments. Preliminary results, using a modified experimental design to improve the time resolution, suggest that the rate of tyrosine phosphorylation can be resolved on a time scale of seconds, with receptor and target proteins reaching peak phosphorylation levels at different time points between 5 and 15 s (data not shown). After 15-30 s, the concentrations of phosphorylated EGFR, Grb2 coprecipitated with EGFR, and phosphorylated PLC␥ started to decrease over time, approaching sustained levels that were much lower than the peaks (Figs. 2 and 3). By contrast, the concentrations of total phosphorylated Shc and Grb2 coprecipitating with Shc increased almost monotonically, reaching high sustained levels (Fig. 3).
We have derived a kinetic model of EGFR signaling pathway to account for the complex cellular responses to EGF. The model has been written in molecular terms as a cascade of protein interactions that involve EGFR, Shc, PLC␥, Grb2, and SOS proteins and phosphorylation/dephosphorylation reactions. Each molecular step of the model is a relatively simple biochemical or physical process. The kinetic parameters have been selected using current information from the extensive literature and/or derived from basic physical-chemical quantities. Our approach provides a quantitative and integrative description of EGFR signal transduction. While the present model simplifies the complexity of the EGFR phosphorylation, it successfully explains the experimentally observed time course of early events in EGF-initiated signaling in liver cells. A transient pattern of the EGFR phosphorylation appears to result from the protection of phosphotyrosine residues against dephosphorylation, as long as a target protein is bound to EGFR. The response behavior is analogous to the propagation of phosphorylation waves through receptor cycles (see Fig. 1). Therefore, the concentrations of target proteins bound to the receptor exhibit transient responses (such as complexes of EGFR with Shc, Grb2, and SOS; see Fig. 5, B and C, lines 1 and 2), whereas the concentration of the phosphorylated Shc and its complexes with Grb and SOS dissociated from the receptor increase almost monotonically (Fig. 5, A and C, lines 3 and 4).
The analysis of sensitivity showed that the response of EGFR signaling pathways to EGF stimulation is stable with respect to changes in the kinetic parameters over a wide range of values. On the other hand, our analysis predicts which factors control the activation state of EGFR and its target proteins. The computational model identifies the kinetic properties of individual reactions that are important for the transient behavior observed experimentally. In particular, computations indicate that after tyrosine phosphorylation of a target protein (Shc, PLC␥) bound to EGFR, its K d for EGFR should increase by at least 1 order of magnitude (compared with the K d for the same unphosphorylated protein). The model demonstrates that the time course of EGFR tyrosine phosphorylation strongly depends on the rate constants of binding and dissociation of phosphorylated Shc protein from EGFR (step 15 in Fig. 1). For instance, a marked decrease in both on and off constants, which leaves the K d of step 15 unchanged, retards the decrease in the total phosphorylated EGFR from the peak to the sustained level. Another prediction of the kinetic modeling is that the phosphatase(s) of EGFR (step 4) must be strongly regulated by the substrate concentration (RP) and cannot be at saturation when transients with pronounced maxima are observed. In other words, the affinity of the phosphatase for the phosphorylated EGFR cannot be very high, so that the range of variation of the concentration of phosphorylated EGFR does not exceed the apparent K m of the phosphatase(s). A surprising prediction is that the effective activities (V max for the saturation condition or V max /K m ratios for a subsaturating (nearly linear) range) should be substantially lower for the phosphatases of phosphorylated Shc and PLC␥ (steps 16 and 8) than for the receptor phosphatase(s) (step 4). A sustained increase in the activities of Shc and PLC␥ phosphatases to the level of the EGFR phosphatase activity will eliminate sharp peaks in the time course of cellular phosphorylation responses to EGF stimulation.
The control of response patterns by the abundance of signaling proteins shown by the model implies that the range of kinetic constants compatible with experimentally observed behavior may change if the assumptions about the amounts of signaling proteins would have to be adjusted on the basis of future measurements. Various cell types and cells in different functional states may show considerable variation in the abundance of signaling proteins. These differences can alter the response patterns to growth factors in a significant manner (cf. Fig. 7).
Importantly, testing the computational results against experimentally observed response patterns restricted model choices and justified some simplifications we made. For instance, the model simplifies EGF binding to EGFR to that described by a single K d . By contrast, some EGF binding studies have been interpreted in terms of two subclasses of receptors with different affinities to EGF (see, e.g. (36,80)). Evidence for two binding sites for EGF was based on the nonlinearity of Scatchard plots, showing negative cooperativity of EGF binding. However, there is no direct experimental evidence for the existence of two stable populations of EGFR with different affinities for EGF. Although several models accounting for nonlinearity of Scatchard plots have been suggested, the precise molecular basis for the difference in affinity of EGF receptors is still unclear. For example, it has been proposed that the high affinity state represents the dimeric form of the receptor (43,(81)(82)(83). However, this model would predict that equilibrium binding data show positive cooperativity, clearly incon-sistent with the experimental observations, which show only neutral or negative cooperativity. In order to account for the experimental data, allosteric binding of EGF (84,85) or the formation of a ternary complex has been proposed where EGF interacts with one or more cell surface molecules that increases binding affinity (86). Indeed, evidence has been found that the high-affinity class of EGF receptors is associated with the cytoskeleton (87,88), specifically with filamentous actin (89). However, an EGFR mutant that lacks the actin-binding site still expressed both high and low affinities for EGF (90). A domain of the intracellular part of the receptor, located within the tyrosine kinase domain, appeared to regulate the affinity for EGF (90). Therefore, an affinity-modulating protein may be a substrate of the EGF receptor kinase. Recently, high affinity EGF binding was found to be lost in HeLa cells overexpressing a mutant dynamin (K44A) (91), suggesting that this affinitymodulating protein interacts with dynamin.
These unresolved issues concerning the different affinities of EGF binding clearly show that future experiments are needed to establish a precise molecular mechanism responsible for different affinity forms of EGFR. However, this knowledge is not of primary importance for modeling the kinetic behavior of EGF-induced signaling, unless the extracellular EGF concentration is well below the K d of high affinity binding, i.e. well below the lowest EGF concentration of 0.2 nM used in our experiments. The K d of 0.6 nM used in the model represent the "pure" receptor lacking the post-binding events. The kinetic scheme (Fig. 1) assumes that all tyrosine-phosphorylated receptor dimers do not release EGF, i.e. the corresponding dissociation constants for EGF became very low (26). Therefore, in this system, the apparent K d for EGF that would be detected in binding experiments should be significantly below the value of 0.6 nM of the receptor per se. Even if the K d for EGF is higher for some of the dimer receptor forms (as in the model proposed in Ref. 85), because of the negative cooperativity of EGF binding only one of the two EGF molecules will dissociate from the receptor. The remaining EGF molecule must have much greater affinity for EGFR (85), so that the amount of free receptor dimer corresponding to any particular phosphorylated EGFR form can be disregarded at EGF concentrations of 0.2 nM or above. Noteworthy, an analysis of the crystal structures of tyrosine kinase receptors suggests that after dimerization and autophosphorylation, an increased kinase activity of the receptor is maintained even when a hormone molecule is removed from the activated receptor dimer (92). Therefore, these receptor dimers will continue to signal, which further supports the assumption of the kinetic scheme in Fig. 1.
Another simplification is that our model considers EGF receptor dimers as a single molecular class. Although it is known that EGFR can form heterodimers with other members of the ErbB growth factor receptor family and initiate diverse downstream pathways (93)(94)(95), considerations of these complex interactions is not required to account for experimental data shown in Figs. 2 and 3. At the same time, our approach provides a strong basis for incorporating distinct subclasses of EGFR in a more complex kinetic model and evaluating their impact on the cellular responses to EGF. The kinetic model suggests which features of the experimental data are indicative of unexpected complexities in the cell's response and also provides a tool to assess the feasibility of different mechanisms to explain anomalous behavior. For instance, the transient phosphorylation response of EGF-stimulated PLC␥ reported in our previous studies (25) is here demonstrated to be indicative of additional regulation of this target protein. The purpose of computer modeling is to provide a basis for guiding experimental analysis and testing explicit hypotheses. A model by itself is not an objective "truth," but it can be used to falsify a specific hypothesis. Therefore, good modeling practice requires a systematic exploration of parameter variance ranges compatible with experimentally observed behavior (96).
The kinetic scheme (Fig. 1) was designed to incorporate only the interactions of EGFR with Shc, Grb2, PLC␥, and SOS (through Grb2). There are other signaling proteins that bind to EGFR in hepatocytes (e.g. phosphatidylinositol 3-kinase, GTPase-activating protein, Eps15), and other, as yet uncharacterized, target proteins may exist. The incorporation of these proteins will generate additional cycles/branches in the kinetic scheme emanating from the phosphorylated receptor-dimer RP (see Fig. 1). Provided the fraction of the receptor bound to additional target proteins is less than the maximal fraction bound to the proteins presented in Fig. 1, these additional interactions will not change the time course of the responses considered here. Moreover, our results make it possible to predict and interpret the time course of EGFR-activated proteins and protein complexes, even those that were not considered explicitly in this paper. In general, transient peaks of the concentrations of complexes of EGFR with activated signaling proteins can be expected (provided rate constants are appropriate). By contrast, the concentrations of activated signaling proteins after their dissociation from the EGF receptor will increase nearly monotonically, unless other inactivation/compartmentation processes occur, such as inhibitory phosphorylation or binding to specific cell structures (with possible subsequent inactivation/degradation). Hence, the kinetics of response patterns of specific target proteins provides information on the nature of downstream regulatory events that can be interpreted in the context of the kinetic scheme presented here.