This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Pribyl, M. Articles by Shvartsman, S. Y. Search for Related Content PubMed PubMed Citation Articles by Pribyl, M. Articles by Shvartsman, S. Y.

Long-Range Signal Transmission in Autocrine Relays

Michal Pibyl^*, Cyrill B. Muratov and Stanislav Y. Shvartsman^*

^* Department of Chemical Engineering and Lewis-Sigler Institute for Integrative Genomics, Princeton University Princeton, New Jersey; and Department of Mathematical Sciences and Center for Applied Mathematics and Statistics, New Jersey Institute of Technology, Newark, New Jersey

Correspondence: Address reprint requests to Stanislav Y. Shvartsman, Dept. of Chemical Engineering and Lewis-Sigler Institute for Integrative Genomics, Princeton University, Princeton, NJ 08544. Tel.: 609-258-4694; Fax: 609-258-0211; E-mail: stas{at}princeton.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS AND DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
Intracellular signaling induced by peptide growth factors can stimulate secretion of these molecules into the extracellular medium. In autocrine and paracrine networks, this can establish a positive feedback loop between ligand binding and ligand release. When coupled to intercellular communication by autocrine ligands, this positive feedback can generate constant-speed traveling waves. To demonstrate that, we propose a mechanistic model of autocrine relay systems. The model is relevant to the physiology of epithelial layers and to a number of in vitro experimental formats. Using asymptotic and numerical tools, we find that traveling waves in autocrine relays exist and have a number of unusual properties, such as an optimal ligand binding strength necessary for the maximal speed of propagation. We compare our results to recent observations of autocrine and paracrine systems and discuss the steps toward experimental tests of our predictions.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS AND DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
The versatility of cell-cell communication relies on sophisticated modules for signal generation, transmission, detection, and processing (Hunter, 2000; Jordan et al., 2000). Compared to signal detection and intracellular transduction, signal generation and transmission are relatively poorly understood. For instance, even though the Epidermal Growth Factor Receptor (EGFR) signaling network has been studied over the past four decades, the molecules that mediate the release of EGFR ligands are being identified only now (Schlessinger, 2000; Lee et al., 2001; Merlos-Suarez et al., 2001; Sunnarborg et al., 2002; Yan et al., 2002). Current work on regulated ligand release focuses on the identification of the relevant molecules and intracellular events. At the same time, it is important to understand how these processes operate in tissues.

Many growth factors and cytokines are released by intracellular or cell surface proteases (Dello and Rovida, 2002). In a number of developmental and physiological contexts, receptors are in excess and their activation is regulated by ligand release and delivery (Freeman and Gurdon, 2002; Kheradmand and Werb, 2002). The availability and activity of ligand-releasing enzymes can be regulated both intra- and extracellularly. Notably, the ligand-releasing enzymes can be activated by the signaling pathways that are stimulated by the ligands they release. Therefore, a positive feedback loop can be established in the regime when the producing cell effectively recaptures and responds to the secreted ligand (Carpenter, 1999; Gschwind et al., 2001; Dello and Rovida, 2002). Such feedback from ligand binding to ligand release is a frequent component of autocrine and paracrine EGFR networks (Dent et al., 1999; Fan and Derynck, 1999; Gechtman et al., 1999; Freeman, 2000; Chen et al., 2001).

Intracellular signaling by the Ras-MAPK pathway can mediate positive feedback in the EGFR system, Fig. 1. In Drosophila, MAPK activated by EGFR induces the transcription of Rhomboid, an intracellular protease, that processes Spitz, a secreted EGFR ligand (Mantrova and Hsu, 1998; Sapir et al., 1998; Wasserman and Freeman, 1998; Guichard et al., 1999; Bang and Kintner, 2000; Hsu et al., 2001; Lee et al., 2001), Fig. 1 A. In mammalian EGFR systems, MAPK can activate the cell-surface ligand-releasing proteases, such as TACE and Kuzbanian (Fan and Derynck, 1999; Merlos-Suarez and Arribas, 1999; Doedens and Black, 2000; Montero et al., 2000, 2002; Merlos-Suarez et al., 2001; Diaz-Rodriguez et al., 2002), Fig. 1 B. The MAPK-mediated feedback in the EGFR system is important in a number of developmental and clinical contexts. The EGFR/MAPK/Rhomboid/Spitz feedback functions throughout the fruit fly development (Casci and Freeman, 1999). The EGFR/MAPK/TACE/TGF network can protect cells against ionizing radiation and prevent the success of cancer radiotherapy (Hagan et al., 2000; Harari and Huang, 2001). These feedbacks rely on the regulation of production and/or activity of ligand-releasing enzymes; other mechanisms relying on the induction of receptor or its ligands have also been described (Doraiswamy et al., 2000; Albanell et al., 2001).

View larger version (31K): [in this window] [in a new window]   FIGURE 1  Positive feedback in autocrine systems. (A) Positive feedback relying on transcription. Intracellular signaling activated by autocrine ligands induces and maintains the transcription of the ligand-releasing enzyme. (B) The feedback is mediated by the activation of the ligand-releasing enzyme. (C) The geometry of cell communication in the model of an autocrine relay. Ligand is diffusing in extracellular medium above the cellular layer. (D) Receptor-mediated processes in the model: ligand release, extracellular transport, reversible binding, and endocytosis. D, ligand diffusivity; kon, ligand-receptor association constant; koff, complex dissociation constant; and kec, ligand-induced internalization rate constant.

A few words on the considered geometry of cell communication. In vitro, this geometry can be realized within the cell and tissue culture assays, where an epithelial layer or a confluent monolayer of autocrine cells is covered by liquid medium; the experimentalist controls the height of the medium (Mandell et al., 2001). This arrangement is also realized in a number of developmental contexts. For example, in Drosophila egg development, a layer of epithelial follicle cells covers a large oocyte (Spradling, 1993; Dobens and Raftery, 2000). Both the oocyte and the follicle cells secrete the EGFR ligands, but EGFR is expressed only in the follicle cells. Hence, ligand diffuses in the gap between the reflective and the receptor-covered surfaces; the size of this gap is less than one micron.

The paper is organized as follows. In the next section we present a model for ligand release, transport, and signaling. We analyze the model using a combination of analytical and computational tools. First, we solve the simplified version of the model in the regime when receptors are in excess and ligand release obeys simple threshold-like kinetics. These assumptions are relaxed during the computational analysis of the model. We conclude with the outline of future directions for modeling and suggest experimental tests of our predictions.

	MODEL FORMULATION

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS AND DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
Equations and scaling
We consider an infinite layer of cells covered by an extracellular medium of height h. The model describes the coupled dynamics of secreted ligands (S), surface receptors (R), ligand-receptor complexes (C), and ligand-releasing proteases (P). The ligand reversibly binds to cell surface receptors; the ligand-receptor complex is degraded through the combination of complex dissociation and receptor-mediated endocytosis. Let us define the ligand concentration at the cellular layer as . Q_r is the rate of receptor synthesis; g_p and g_r denote the linear gains in the protease and ligand production, respectively. The balance equations and the corresponding boundary conditions take the following form:

(1)

(2)

(3)

(4)

(5)

where t is time, x is the spatial coordinate parallel to the cellular layer, and y is the spatial coordinate perpendicular to it (see Fig. 1 C).

The boundary conditions at the surface of the cellular layer account for diffusive flux of extracellular ligand, reversible ligand-receptor binding, and generation by intracellular or cell surface ligand-releasing proteases. In our model, ligand generation is controlled by the activity of the protease and has zero order with respect to the amount of cell-associated ligand precursor. In other words, the cell-associated ligand precursor is in excess and is not significantly depleted by the proteolytic release. This is true in at least two experimental EGFR systems: Spitz release by Rhomboid in Drosophila development (Rutledge et al., 1992; Freeman, 1997; Stemerdink and Jacobs, 1997; Stevens, 1998; Wasserman and Freeman, 1998; Bang and Kintner, 2000) and release of TGF from autocrine epidermoid carcinoma cells (Dent et al., 1999).

In the present formulation of the model, ligand-receptor complex is internalized in the first-order process with rate constant k_ec. A straightforward extension of the model can account for the following processes of ligand/receptor trafficking. This will require introduction of additional state variables corresponding to species at various stages of endocytic cycle (Burke et al., 2001; Waterman and Yarden, 2001). The parameters of the model and their typical values are listed in Table 1.

We assume that protease synthesis (or activation) is an algebraic function of the number of ligand-receptor complexes. The nonlinearity (C) combines a large number of processes leading from receptor binding to protease generation. This simple model assumes that these steps are fast, and that protease generation can instantaneously follow the number of occupied receptors. At this stage, we need this assumption to advance our analysis. Toward the end of the paper, we discuss and computationally analyze the potential role of a time lag between receptor occupancy and receptor activation (Fig. 8 B). In the simplest case, this lag can be parameterized by a constant time delay. This would modify the right-hand side of Eq. 4 (the sigmoidal function would depend on the number of complexes at a previous moment in time: (C(t - t_D)). In terms of the standard control theory, the model for protease dynamics is then classified as a first-order system with delay (Ogunnaike and Ray, 1994).

View larger version (12K): [in this window] [in a new window]   FIGURE 8  Robustness of traveling fronts with respect to several modifications in the model. The relationship between the dimensional (V) and the dimensionless (v) velocity is given by V = vLxkp. (A) Dimensionless velocity monotonically decreases with growing value of the Hill coefficient n in the sigmoidal nonlinearity. For n = 2, the difference between dimensionless velocity of the full model and the velocity of the simplified model with sharp nonlinearity of protease activation in the ligand-limited regime (n +, c 0, 0) is 25%. For n = 5, the difference is only 4%. (B) Velocity of wave propagation monotonically decreases with increasing time delay in the sigmoidal nonlinearity of protease activation. Dimensionless time delay D is plotted with units 1/kp (when kp = 0.01 min-1). The figure means that the wave propagation will not stop for the time delay 500 min and less. (C) Increase of the relative time scale of ligand-receptor complex processes c slows down the front. The other dimensionless model parameters and parameters of Hill nonlinearity in the studies (A) and (C) were kept at the same values as in Fig. 6, A and B; in the study, (B) used the same values as in Fig. 5 A, and = 0.1.

(6)

where

(7)

The time scale is set by the dynamics of protease degradation. The vertical length scale is scaled by the height of the medium. The horizontal length scale is derived from the previous analysis of the spatial extent of autocrine loops. Specifically, the length scale L_x appears in the expressions for the cumulative distribution functions that characterize the extrema of random trajectories followed by ligands diffusing above the receptor-covered plane (Shvartsman et al., 2001).

After rescaling, we get:

(8)

(9)

(10)

(11)

(12)

There are six dimensionless parameters in the rescaled model:

(13)

The first of these, a, is the Damköhler number that characterizes the relative rates of diffusion and binding (Deen, 1998). Note that a is also the ratio of the geometric and dynamic length scales, h and L_x, respectively. The second dimensionless group, ß, characterizes the relative rates of endocytosis and dissociation; is proportional to the ratio of the ligand and receptor generation rates (i.e., ratio of g_rg_p/k_p and Q_r); _s, _r, and _c characterize the relative time scales of the extracellular ligand, and of the free and bound receptors, respectively.

Steady-state multiplicity and traveling fronts
The spatially uniform steady states of Eqs. 8–12 satisfy:

(14)

The first equation is the balance of linear and sigmoidal functions, a classical mechanism for the generation of steady state multiplicity; Fig. 2 A (Keener, 1998; Fall et al., 2002). Thus, depending on the parameters of the sigmoidal nonlinearity, the system of Eq. 14 may have either one or two stable steady states. In the case of bistability, the first of the steady states corresponds to the "off" state of the autocrine loop, with the zero rate of ligand release—c = 0, p = 0, r = 1, s = 0. The second steady state corresponds to the "on" state of the autocrine loop, with nonzero level of ligand-receptor complexes and appreciable rate of ligand release.

View larger version (13K): [in this window] [in a new window]   FIGURE 2  (A) Graphical solution of Equation 14, where (c) is chosen to be a Hill function (c) = cn/(kn + cn). (B) The case of the Heaviside nonlinearity: (c) = H(c - c0). No ligand is released in the "off" steady state. Autocrine loop is activated in the "on" steady state.

Our main finding is that such fronts may indeed be realized in autocrine relay systems. They arise as a result of the coupling of cells by secreted autocrine ligands; this coupling is provided through the diffusion of the ligand in the extracellular space. Note that this form of coupling is different from the one that is usually encountered in models of bistable media, in which an autocatalytic substance both propagates by diffusion and undergoes the autocatalytic reaction at the same time (Mikhailov, 1994; Keener, 1998). In our model, the diffusing variable is just a messenger mediating communication within the layer of autocrine cells. The positive feedback results from the interaction between the messenger and the localized species—in our case, bound receptors and ligand-releasing proteases. Furthermore, the nature of transport in our problem makes the coupling nonlocal. A localized change in the values of nondiffusing variables affects the entire distribution of the messenger concentration.

The fronts are analyzed using a combination of analytical and computational approaches. In the computational experiments, we were starting with the autocrine layer in the "off" state and following the spatiotemporal evolution of the localized disturbance that activated the protease dynamics. We characterize the shape and speed of the traveling front that is established after the initial transient period. In the following, the description of computational experiments is preceded by the discussion of the analytical solution for the fronts in a certain limiting regime. Specifically, we consider the regime that is characterized by the fast dynamics of ligand-receptor complexes, the vast excess of empty receptors, and infinitely sharp nonlinearity of the protease production term. The analytical solution obtained in this regime enables a complete parametric analysis of the fronts and reveals a number of their unusual properties.

Analytical solution of the simplified model
The ligand-limited regime is characterized by the small ratio of the rates of ligand release and receptor synthesis, << 1. In this case, the first term in the right side of Eq. 9 for the dynamics of empty receptors can be neglected. As a result, the balance for empty receptors is effectively uncoupled from other variables and the dimensionless concentration of free receptors is constant: r 1. The magnitude of characterizes the extent to which receptor is in excess with respect to ligand (see the discussion of the dimensionless groups in the previous section), and is not explicitly related to the time scales of binding. In the large number of developmental and physiological settings, receptor is in excess and its activation is controlled by ligand availability (Freeman and Gurdon, 2002; Kheradmand and Werb, 2002; Sunnarborg et al., 2002). When receptors are in vast excess, their level is not strongly affected by the ligand-induced endocytosis and can be assumed constant across the traveling front.

When the time scale characterizing the equilibration of ligand-receptor complexes in the surface layer ( = 0) is fast, and _c << 1, these species can be eliminated using a pseudo-steady-state approximation: c(,) = s(,0,). The time scale characterizing the relaxation of ligand-receptor complexes is expected to be short compared to the time scale characterizing the protease activation. This time scale is set by the magnitude of the rate constant k_p, which determines the rate with which the protease reacts to a change in the level of its activation. For the case of Rhomboid, an intracellular ligand-releasing protease that is transcriptionally activated by the EGFR signaling (Mantrova and Hsu, 1998; Wasserman and Freeman, 1998; Urban et al., 2002), this time scale of the protease activation is definitely longer than the one for binding. In the case of TACE, a cell surface metalloprotease activated by the processes that does not rely on transcription, this assumption has to be checked experimentally (Dent et al., 1999; Gechtman et al., 1999; Dello and Rovida, 2002).

Thus, under the conditions _c >> 1, << 1, the system 8–12 is reduced to the two equations for the extracellular ligand (s), Eq. 8, and the ligand-releasing protease (p), Eq. A1, in the Appendix. The only source of nonlinearity is provided the sigmoidal generation term in protease dynamics. A sharp sigmoidal nonlinearity can be approximated by a Heaviside function: (c) H(c - c₀), Fig. 2 B. In this case, the problem can be solved analytically; see Appendix. The main result is the implicit equation that links the front speed to the geometric, molecular, and cellular parameters of the autocrine loop, Eq. A18. This equation was solved graphically to produce the parametric plots in the next section.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS AND DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
Fronts in the asymptotic regime
In the ligand-limited regime ( << 1) with fast binding (_c << 1) and sharp nonlinearity ((c) H(c - c₀)), the properties of propagating fronts depend on four dimensionless parameters. The value of c₀ defines the threshold of the nonlinearity, _s is the relative time scale for the dynamics of the extracellular ligand, compares the rates of binding and transport, and ß depends on the kinetic parameters of the ligand-receptor complex. Given the values of these parameters, the formulas derived in the Appendix provide the spatial profile and the speed of the propagating front. Fig. 3 presents the spatial profiles of ligand and protease for a particular set of parameters corresponding to the advancing "on" state of the autocrine loop.

View larger version (11K): [in this window] [in a new window]   FIGURE 3  A traveling front in the ligand-limited regime with fast binding and sharp nonlinearity of protease activation (see Appendix A). (A) and (B) present the profiles of the extracellular ligand and ligand-releasing protease, respectively. The arrows denote the direction of wave's motion. The value of the dimensionless velocity is v = 1.967; the threshold in the Heaviside function, c0 = 0.25. Other dimensionless model parameters are = 55.36, ß = 0.5, s = 5.439 x 10-4 (based on kon = 1 x 102 µM-1 min-1, R0 is based on 5 x 104 receptors/cell surface area (25 µm2), kec = 0.1 min-1, koff = 0.1 min-1, D = 1 x 10-10 m2 s-1, kp = 0.01 min-1, h = 1 x 103 µm).

Figs. 4 and 5 show the effects of molecular, cellular, and geometric parameters on the speed of the propagating solutions. We fix the diffusion coefficient and the height of the extracellular medium and analyze the effect of the binding and trafficking parameters. The height of the medium is 1 mm, a typical value encountered in cell and tissue culture assays of autocrine loops (DeWitt et al., 2001); ligand diffusivity is set to 10^-6 cm²s^-1, a typical value for the diffusivity of a protein in solution.

View larger version (40K): [in this window] [in a new window]   FIGURE 4  Dependence of the dimensional front speed on the kinetic and transport parameters in the ligand-limited regime with fast binding and sharp nonlinearity of protease activation (see Appendix A). (A) Dependence on the rate of forward binding on the forward-binding rate constant, konR0, ligand-induced internalization rate constant kec and complex dissociation constant koff. The dark red area corresponds to the maximal value of the velocity (v 4 x 10-2 µm s-1). Conversely, the dark blue area denotes the minimal value of the velocity (V 4 x 10-4 µm s-1). The threshold in the Heaviside function c0 = 0.25. Other model parameters are h = 1 x 103 µm, D = 1 x 10-10 m2 s-1, kp = 0.01 min-1. (B) A one-dimensional cut (kec/(kec + koff) = 0.5) through the two-parameter plot in (A). The arrows indicate that increasing the ratio kec/(kec + koff) shifts the velocity maximum to a lower value of konR0.

View larger version (32K): [in this window] [in a new window]   FIGURE 5  Dependence of the wave propagation velocity on the ratio of the geometric length scale h and the dynamic length scale Lx in the ligand-limited regime with fast binding and sharp nonlinearity of protease activation (see Appendix A). The relationship between the dimensional (V) and the dimensionless (v) velocity is given by V = vLxkp. (A) Dependence of the dimensionless velocity on the ratio h/Lx computed for several thresholds c0 in the Heaviside function. The propagation velocity increases with the ratio until log10(h/Lx) >> 1. Above this value, the velocity does not depend on this geometric parameter. The parameters are ß = 0.5, s = 5.439 x 10-4 (based on kon = 1 x 102 µM-1 min-1, R0 is based on 5 x 104 receptors/cell surface area (25 µm2), kec = 0.1 min-1, koff = 0.1 min-1, D = 1 x 10-10 m2 s-1, kp = 0.01 min-1). (B) Color maps of the distributions s(, ) of the ligand in a traveling wave in both horizontal and vertical directions, with c0 = 0.25. The spatial distributions of the ligand correspond to the filled circles in (A). The profiles were computed for log10(h/Lx) = -1, 0, 1, and 2. The dark/bright regions correspond to the on/off steady states (s = 1, s = 0), respectively.

The nonmonotonic dependence of the front speed on the number of surface receptors and/or the forward binding rate constant is a result of the particular nature of coupling in autocrine relays. A sufficient number of cell surface receptors must be occupied to sustain the protease activation; this explains the ascending part of the curve. At the same time, the traveling front is slowed down when ligand spends most of its time bound to cell surface receptors; this qualitatively explains the descending part of the curve. A more detailed study of the parametric dependences of the velocity will be presented elsewhere. The maximum in the dependence of the front speed on k_onR₀ is the robust property of autocrine relays; it exists for any value of the threshold in the nonlinearity and the rate constants of dissociation and internalization. The speed of the front depends on the kinetic properties of ligand-receptor binding. Hence, it is not possible to compare front speeds in ligand-receptor systems based only on the corresponding equilibrium binding constant.

The height of the medium can be easily adjusted in cell and tissue culture assays. Fig. 5 A presents the dependence of the front speed on the ratio of the geometric and dynamic length scales in the problem (i.e., ratio of h and L_x). When the ratio of these two length scales is small, the front speed is an increasing function of the medium height. Once again, this is the result of the nature of coupling in autocrine epithelial layers. For small height of the medium, the time interval between successive binding events of secreted ligand is short; the ligand spends a large fraction of its time in the bound form. Consequently, increasing the height of the medium increases both the duration and the span of the ligand trajectories and, as a result, the speed of the propagating front. In this regime, the fronts are essentially one-dimensional, in a sense that they are nearly uniform in the vertical direction; see the two top profiles in Fig. 5 B. This is easy to understand: the small value of the medium height translates into the low value of the Damköhler number, , that determines variation of the solution in the vertical direction; see Eqs. 8–12. In fact, a conventional "thin fin" approximation (Deen 1998) derived in the limit 0 captures these fronts and their dependence on the height of the medium very accurately (results not shown).

When the height of the medium is comparable to the dynamical length scale L_x of the problem, the fronts become fully two-dimensional; see the two bottom profiles in Fig. 5 B. The dependence of the front speed on the medium's height h asymptotically approaches a constant value, Fig. 5 A. This behavior can be explained in terms of the results of a recent analysis of autocrine loops in the semiinfinite medium (Shvartsman et al., 2001). The statistical properties of both the vertical and the lateral spans of the random trajectories followed by autocrine ligands are characterized by the single length scale given by L_x. Specifically, half of the ligand trajectories are captured for the first time before their maximal lateral displacement reaches L_x. This is precisely the dynamical length scale in our problem; see Eq. 6. In other words, autocrine loops are spatially localized even in the semiinfinite medium. This means that, above some critical value, increasing the height of the medium would stop contributing to the effective range of the messenger variable in our problem. Consequently, the dependence of the front speed on the height of the medium should vanish at high ratios of the geometric and the dynamic length scales in the problem, Fig. 5 A.

Numerical simulations of propagating fronts
To test the predictions of the asymptotic analysis, we carried out numerical simulations of the full model described by Eqs. 8– 12. The spatial derivatives in Eqs. 8 and 12 were replaced by finite differences and the resulting set of ODEs was solved using an efficient explicit solver for parabolic partial differential equations (Sommeijer et al., 1998).

Our main result can be summarized as follows: traveling fronts predicted in the simplified model exist for a wide range of parameters in the full model. Moreover, they exhibit the properties predicted by the analysis of the simplified model.

Fig. 6, A and B presents the profile of the fully developed front computed in the full model with the Hill nonlinearity, a nonzero value of the time scale of ligand-receptor complex, and a nonnegligible ratio of the ligand and receptor generation rates. This front was computed in a transient simulation that started with the system in the "off state" and followed the evolution of a localized disturbance activating the ligand release. After the initial transient, this disturbance evolves into a fully developed self-similar solution—a constant speed front, Fig. 6 C. Fig. 6 D shows that, in agreement with the predictions of the analysis of the ligand-limited regime with fast binding and sharp nonlinearity of protease activation, the speed of the fronts in the full model exhibits a maximum with respect to the intensity of the binding process (k_onR₀).

View larger version (21K): [in this window] [in a new window]   FIGURE 6  Numerical analysis of the full model. (A) and (B) Ligand and protease distributions computed across the front. The ligand release function is modeled by a Hill function (c) = cn/(kn + cn), with k = 0.25 and n = 5. Other parameters are = 1, ß = 0.5, = 0.2, c = 0.1, s = 3.334, and r = 2 (based on kon = 1.807 µM-1 min-1, R0 = Qr/ker is based on 5 x 104 receptors/cell surface area (Qr = 500 receptors/cell surface/min, ker = 0.01 min-1, cell surface area 25 µm2), kec = 0.1 min-1, koff = 0.1 min-1, D = 1 x 10-10 m2 s-1, kp = 0.02 min-1, h = 1 x 103 µm, Ql = grgp/kp = 50 ligand molecules/cell surface/min). (C) The wave velocity was computed from the linear part of the dependence of the front position on time. (D) Dependence of the wave propagation velocity on the ligand-receptor affinity konR0. The points are the results of the numerical analysis; the curve is drawn to guide the eye. The values of the parameters of the Hill nonlinearity and the other dimensional model parameters (except kon and R0) are the same as in (A) and (B).

Fig. 7 illustrates that traveling waves in the model of autocrine relays can "jump" over the heterogeneities in the cellular layer. This behavior is not unexpected, inasmuch as the mechanism of signal transmission is mediated by diffusion through the extracellular medium. What is interesting, however, is that the fronts can be both accelerated and impeded by the heterogeneities. This effect can be explained by the maximum in the dependence of the front speed on the number of cell surface receptors, Fig. 6 D. When the parameters of the monolayer correspond to the ascending part of the curve, the front is slowed down by the heterogeneity. In this regime, the heterogeneity provides a localized break in the positive feedback loop, Fig. 7 B. However, to the right of the maximum, when the front speed is a decreasing function of the number of receptors, the front will be accelerated by the local heterogeneity. In this regime, the front motion in the unperturbed layer is slowed down by frequent ligand-receptor binding. Thus, heterogeneity characterized by the absence of cell surface receptors provides a localized "escape path" for autocrine ligands, increasing the range of their action, Fig. 7 D. Notice that in both cases the fronts regain the constant velocity after their interaction with local heterogeneity. We have found that for the front to "jump across the bump," the size of heterogeneity has to be of the order of the dynamic length scale in the problem, L_x = (Dk_er)/k_onQ_r), Fig. 7 C. These computational experiments underscore the unusual properties of traveling waves in autocrine systems. We can also predict the effect of other heterogeneities, e.g., in the value of ligand diffusivity. For example, in the regime when ligand binding slows down the front (the descending part of the curve in Fig. 4 B), a "defect" with increased diffusivity will speed up the front.

View larger version (26K): [in this window] [in a new window]   FIGURE 7  Interaction of fronts with heterogeneities in the cellular layer. (A) The heterogeneities can both accelerate and slow down the fronts. The parameter x shows an extent of the heterogeneity. (B and D) The location of the front as a function of time in an intact (solid line) and damaged layer (broken line). (B) The case of low binding rate, konR0 = 0.1 µms-1 (based on kon = 1.807 µM-1min-1). The extent of the damaged tissue is x/Lx = 2.5. The values of the parameters of the Hill nonlinearity and the other model parameters are the same as in Fig. 6, A andB. (C) An increase of the heterogeneity extent stops the traveling front. It happens for x/Lx 6.5 in the case of the low binding rate; see (B). (D) The case of high ligand-receptor affinity konR0 = 10 µms-1 (based on kon = 1.807 x 102 µM-1min-1). The extent of the damaged tissue is x/Lx = 10. Other parameters are as in Fig. 6, A and B.

The _c << 1 assumption can be readily relaxed. Fig. 8 C presents the dependence of the front speed on this parameter. The general effect of increasing the time scale for ligand dynamics is to slow the front speed down.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS AND DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
We have developed a mechanistic model of autocrine relays. The model enables a systematic analysis of nonlinear interaction between ligand release, transport, and binding. Our transport model accounts for the nonlocal coupling between different regions of the epithelial layer; the range of the coupling is directly related to the transport, binding, and endocytosis of the secreted ligand. When combined with a positive feedback by ligand release, the model robustly predicts the existence of propagating fronts. These fronts travel at constant speeds and can deliver signals to arbitrary distances. This wave-like propagation is very different from the one that relies on pure diffusion. Specifically, we predict the existence of traveling waves mediated by the autocrine release of growth factors. In other systems, traveling waves can be mediated by the regulated release of extracellular calcium, ATP, and cAMP (Kessler and Levine, 1993; Tang and Othmer, 1995; Tang et al., 1996; Palsson et al., 1997; D'andrea et al., 1998; Muller et al., 1998; Sneyd et al., 1998; Homolya et al., 2000; Scemes et al., 2000; Dormann et al., 2001; Shiga et al., 2001; Hofer et al., 2002; Schuster et al., 2002). The fronts mediated by autocrine loops have a number of unusual properties, such as a nonmonotonic dependence of the front speed on the total number of surface receptors and/or the forward binding rate constant. Wave-like patterns of intercellular communication can rely on multiple mechanisms. In several systems, a combination of paracrine and gap-junctional communication has been described (Hirata et al., 1998; Sauer et al., 2002). The physiological significance of these waves can range from robust and long-range signal delivery in normal tissue physiology (e.g., in response to injury; Mandell et al., 2001; Katz et al., 2002) to wave-mediated patterning of epithelial tissues (a propagating wave can leave a patterned state behind.

The predicted velocities of fronts in paracrine relays are dictated by the time scale of protease activation and by the effective diffusivity determined by the combination of the ligand diffusion coefficient and the cell surface receptor density. The velocity of waves in autocrine relays mediated by secreted growth factors is lower than the velocity observed in axonal propagation and in the intracellular waves. On the other hand, these waves might operate in the processes of developmental pattern formation (Freeman, 2000; Freeman and Gurdon, 2002) and tissue remodeling in response to injury. At this point, one has to be careful in comparing our predictions to the observations of fronts in paracrine relays in growth factor systems. The reported experiments were done in large dishes, where medium convection could contribute to the increase of the apparent propagation velocity (Mandell et al., 2001; Katz et al., 2002).

We have considered only planar fronts. Our preliminary results indicate that convex/concave fronts will move slower/faster than the planar fronts for the same values of kinetic/transport parameters. This is expected from the behavior of conventional trigger waves in reaction-diffusion systems (Mikhailov, 1994). Hence, we expect the planar interfaces to be stable with respect to the lateral perturbations. In the future, it will be interesting to determine if autocrine relays can support fully two-dimensional fronts, such as the cusp-like interfaces described in (Brazhnik and Tyson, 1999). Recent experiments on paracrine waves in truly three-dimensional geometries, such as tumor spheroids, make the multidimensional extension of our analysis particularly relevant (Sauer et al., 2002).

At this point, there is no direct evidence for traveling waves mediated by the regulated secretion of peptide growth factors. Recently, however, a paracrine mechanism has been proposed to cause the wave-like spread of signaling activity induced by localized mechanical stimulus (Mandell et al., 2001). In addition, long-range cell communication mediated by the interaction between paracrine growth factors and ERK2 signaling has been observed in cell culture experiments with wounded osteoblast monolayers (Katz et al., 2002). There, fibroblast growth factors have been proposed as paracrine signals. A wealth of experiments supporting the release-mediated positive feedback in growth factor systems indicates that these waves can be observed (Hunter, 1998; Carpenter, 1999; Dent et al., 1999; Fan and Derynck, 1999; Gechtman et al., 1999; Moghal and Sternberg, 1999; Albanell et al., 2001; Chen et al., 2001; Gschwind et al., 2001; Ma et al., 2001; Merlos-Suarez et al., 2001; Murphy et al., 2001; Carraway and Sweeney, 2002; Diaz-Rodriguez et al., 2002; Montero et al., 2002). Usually, this feedback is studied in cell culture, where cells are randomly dispersed on the surface. We predict that experiments with confluent cell monolayers or with model epithelial layers will uncover traveling waves mediated by autocrine growth factors. The EGFR network is a particularly promising experimental system. The experiments require an ability to locally stimulate ligand release and to follow the progress of activation across the cellular layer. Localized stimulation can be mechanical, chemical, or rely on gamma or UV radiation. Excitation can be followed either in situ (e.g., using fluorescent reporters) or ex situ (e.g., using Western blot epithelial layers at different time points). The A431 carcinoma cell line with a positive EGFR/MAPK/TGF feedback is a good experimental system for these studies (Dent et al., 1999; Shvartsman et al., 2002).

Almost invariably, nonlinear traveling waves in cell communication are modeled as reaction-diffusion systems, where the amplification processes in the active medium happen in the same region as transport. In our model, these processes are separated in a sense that active processes at different parts of the surface are nonlocally coupled by diffusion through the volume. The properties of traveling waves in such systems are yet to be characterized in details. From the mathematical standpoint, nonlinear waves in multiphase reaction-transport systems were studied by Othmer (1975).

Traveling waves in autocrine and paracrine relays can propagate through the discontinuities in the phase that generates the active signal. In vitro, this property can be revealed by experiments with "scratched" monolayers, as was done in an elegant study by Mandell (et al., 2001). In vivo, this setup can be realized by generating the clones of cells lacking the receptor for the paracrine ligand; see Tabata (2001) for the recent review of relevant genetic techniques. In both cases, experiments should be spatially resolving and cannot rely just on localized measurements.

The positive feedback from ligand binding to ligand release is just one of the many regulatory patterns in autocrine and paracrine networks. Some of the other mechanisms involve binding-stimulated receptor shedding (Ni et al., 2001), decrease of receptor mRNA stability (Sturtevant et al., 1994), and release of extracellular ligand/binding components (Guan et al., 2001; Freeman and Gurdon, 2002; Gerlitz and Basler, 2002). The physiological significance of these processes is only starting to be appreciated.

The next stage of modeling should account for the dynamics of intracellular signaling. A number of "lumped" models of EGFR signaling can be incorporated into our models of autocrine epithelial layers (Bhalla and Iyengar, 1999; Kholodenko et al., 1999; Brightman and Fell, 2000; Haugh et al., 2000; Schoeberl et al., 2002). For this, a model of protease dynamics is necessary. A simplified description, accounting for the stimulation-induced downregulation of the protease and the kinetic patterns of EGFR ligand release, has been recently developed (Dong and Wiley, 1999; Fan and Derynck, 1999; Gechtman et al., 1999; Doedens and Black, 2000; Shvartsman et al., 2002). A growing number of studies of regulated ligand release can be used to develop and validate a more detailed model. This model would simultaneously account for multiple species and processes leading to binding-induced ligand release (Albanell et al., 2001; Merlos-Suarez et al., 2001; Umata et al., 2001; Zhang et al., 2001; Diaz-Rodriguez et al., 2002; Kheradmand and Werb, 2002; Montero et al., 2002; Sunnarborg et al., 2002).

We have emphasized the positive feedback and activation processes in autocrine relays. This is sufficient to capture the leading part of the propagating wave. Various intracellular processes, working on the time scales of minutes to hours, may switch the system off. This would convert the traveling front to a traveling pulse. When the deactivation processes are slow, the traveling pulse has a well-defined front edge (Meron, 1992), and pulse speed is well-approximated by the analysis that neglects deactivation in the back of the wave. Increasing the ratio of the time scales of deactivation and activation will eventually destroy the traveling pulse. In the case of conventional reaction-diffusion systems, this transition is through the saddle-node bifurcation (Mikhailov, 1994). Among the processes that contribute to the deactivation in autocrine relays are well-established adaptation in growth-factor-induced ERK2 signaling, downregulation of level of ligand precursor, and downregulation of the protease activity. The last process has been clearly demonstrated in the recent experiment of Doedens and Black (2000). The TACE protease that regulates the availability of EGFR ligands (Sunnarborg et al., 2002) is activated by a number of intracellular pathways. The activated form of this enzyme then cleaves the ectodomains of transmembrane molecules. At the same time, the active form of the enzyme is primed for degradation and is cleared from the surface through the endocytic pathway.

	APPENDIX

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS AND DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
The analytical solution is constructed in the limiting case of _c 0, 0, and (c) = H(c - c₀), where H(x) is the Heaviside function and c₀ is the threshold. Under these assumptions, we have r = 1 and on the time scale of the wave, so Eqs. 9–12 simplify to

(A1)

(A2)

We are looking for the traveling wave solutions s = s( - v,), p = p( - v) connecting the steady state s = 1, p = 1 at = - with the steady state s = 0, p = 0 at = +, moving with speed v > 0 to the right. Introducing the reference frame moving with the wave

(A3)

and rewriting Eqs. 8 and A1 in terms of , we obtain

(A4)

(A5)

where we assumed that s(,0) is a monotonically decreasing function of and placed the point of the onset of the protease production at the origin.

Eq. A5 can be straightforwardly integrated:

(A6)

This solution can be substituted into Eqs. A2 and A4, which can then be solved by separation of variables. Let us look for the solution of Eqs. A2 and A4 in the form

(A7)

where _n() are the orthonormal eigenfunctions of the Sturm-Liouville problem

(A8)

A simple calculation gives

(A9)

(A10)

Substituting Eq. A7 into A4, multiplying it by _n, and then integrating over , we obtain an equation for , where and are the solutions for > 0 and < 0, respectively:

(A11)

This equation has to be supplemented by the boundary conditions at = 0 and ±. At infinity we must use the exact fact that s(+,0) = 0 and s(-,0) = 1. We should also use the fact that s together with its first derivatives must be continuous at = 0. In view of Eq. A7, this translates to

(A12)

(A13)

Once again, multiplying Eqs. A12 and A13 by _n, ntegrating over , and combining them with the behavior at infinity, we obtain that , which should satisfy the following boundary equation,

(A14)

(A15)

After a rather long, but straightforward calculation, we obtain

(A16)

and

(A17)

Substituting this solution into Eq. A7 and using it in the self-consistency condition in Eq. A5, we obtain the implicit expression for the speed of the traveling wave,

(A18)

	ACKNOWLEDGEMENTS

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS AND DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
We thank Giovanni Samaey for help with the DDE-BIFTOOL software and Jim Mandell for many helpful discussions.

Our work has been supported by the National Science Foundation.

Submitted on September 17, 2002; accepted for publication November 12, 2002.

	REFERENCES

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS AND DISCUSSION CONCLUSIONS APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
Albanell, J., J. Codony-Servat, F. Rojo, J. M. Del Campo, S. Sauleda, J. Anido, G. Raspall, J. Giralit, J. Rosello, R. I. Nicholson, J. Mendelsohn, and J. Baselga. 2001. Activated extracellular signal-regulated kinases: association with epidermal growth factor receptor/ transforming growth factor a expression in head and neck squamous carcinoma and inhibition by anti-epidermal growth factor treatments. Cancer Res. 61:6500–6510.[Abstract/Free Full Text]

Bang, A. G., and C. Kintner. 2000. Rhomboid and Star facilitate presentation and processing of the Drosophila TGF-alpha homolog Spitz. Genes Dev. 14:177–186.[Abstract/Free Full Text]

Bhalla, U., and R. Iyengar. 1999. Emergen