| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
ibyl *
* 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, Department 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 INTRODUCTIONDESCRIPTION OF THE MATHEMATICAL...RESULTSAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONDESCRIPTION OF THE MATHEMATICAL...RESULTSAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
). In general, cell-to-cell communication produces spatially nonuniform patterns in the expression of genes that guide the development of tissues and organs. The design principles of epithelial patterning are being formulated only now (Freeman and Gurdon, 2002). An important family of epithelial patterning mechanisms relies on secreted chemical signals. Typically, a ligand released by a group of cells interacts with the extracellular matrix and cell surface receptors as it spreads through the tissue. Ligand transport can be integrated with positive and negative intracellular feedback loops (Freeman, 2000). For example, ligand-receptor binding can stimulate ligand synthesis and secretion (Freeman and Gurdon, 2002). Ligand release can be regulated by the occupancy of cell surface receptors. Receptor occupancy, in turn, may be determined by the balance between ligand transport, binding, and degradation; see Fig. 1 D.
|
), where pattern formation proceeds in the follicular epithelium—a layer of columnar epithelial cells that envelop the oocyte. Follicle cells are much smaller (5–7 µm) than the oocyte (100–300 µm); see Fig. 1 A. Reciprocal oocyte/follicle cell interactions pattern the eggshell and establish the embryonic axes (Van Buskirk and Schupbach, 1999). These events rely on the epidermal growth factor receptor (EGFR), a well-studied receptor tyrosine kinase (Nilson and Schupbach, 1999; Wells, 1999). EGFRs are uniformly distributed across the follicular epithelium and are absent on the oocyte surface (Sapir et al., 1998). EGFR is activated by ligands secreted from the oocyte and from the follicle cells themselves. Secreted ligands diffuse in the thin (<1 µm) gap between the follicle cells and the oocyte; see Fig. 1 B. Release of EGFR ligands is regulated by the intracellular proteases (Urban et al., 2002). Interestingly, the expression of these proteases (Rhomboids), is positively regulated by EGFR signaling (Hsu et al., 2001; Peri et al., 1999; Sapir et al., 1998; Wasserman and Freeman, 1998); see Fig. 1 E. To summarize, ligands diffuse in a thin gap between receptor-covered epithelium and a reflective surface; receptor activation stimulates further ligand release by activating the expression of the intracellular protease.
Similar signaling/transport arrangements are encountered later in fruit fly development and in other species, both for the EGFR and other signaling systems (Doraiswamy et al., 2000; Freeman and Gurdon, 2002).
The gene expression patterns in developing epithelial layers can be very fine-grained. In many cases, the width of the signaling patterns is only a couple of cells (Carmena et al., 2002; Hatini and DiNardo, 2001; Peri et al., 1999; Ruohola-Baker et al., 1993; Sapir et al., 1998). What is the appropriate biophysical description for such systems? There have been several attempts to use continuum models for the analysis of such patterns, e.g., those of Lander and co-workers (Lander et al., 2002), and Shvartsman and co-workers (Shvartsman et al., 2002). These models assume that the relevant length scale of the pattern is greater than the size of a single cell. The validity of such an assumption may be difficult to reconcile with the fine-grained nature of experimentally observed signaling patterns. Starting from the work of Othmer and Scriven, a number of discrete models for cellular layers has been proposed (Collier et al., 1996; Meir et al., 2002; Monk, 1998; Othmer and Scriven, 1971; Owen et al., 1999, 2000; von Dassow et al., 2000). However, in these models the form of cell-to-cell couplings has been chosen on purely phenomenological grounds.
Here, we systematically derive discrete models of cell-to-cell communication from a mechanistic description of autocrine and paracrine signaling in epithelial layers. These models are discrete because they treat each cell individually. The models are also long-ranged: they use the state of the entire layer in describing the dynamics of each cell and thus take into account non-nearest neighbor cell-to-cell interactions. Our derivation is based on the introduction of coupling coefficients that characterize the communication of cells by secreted ligands. These coupling coefficients are directly linked to the biophysical parameters of the transport problem. Central to our approach is the use of the separation of timescales between diffusing and intracellular species. We argue that binding and transport are fast and are therefore dynamically slaved to the slow intracellular variables when the intracellular processes involve transcription and protein synthesis. With the coupling coefficients at hand, we are then able to formulate discrete models that account for particular intracellular processes, such as receptor-mediated activation of ligand release, and study their properties associated with discreteness.
|
DESCRIPTION OF THE MATHEMATICAL MODEL |
|---|
|
TOPABSTRACTINTRODUCTIONDESCRIPTION OF THE MATHEMATICAL...RESULTSAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
). Following the standard receptor binding analysis, this approximation requires that the concentration of ligand is less than the equilibrium binding constant, kD = koff/kon. With this in mind, a mechanistic model of transport and signaling in an idealized epithelium consisting of a two-dimensional periodic array of identical cells in a single flat layer geometry takes the following form (Fig. 1, B and D):
 | (1) |
 | (2) |
 | (3) |
 | (4) |
is the concentration of ligand in the extracellular space, 
is the number of ligand-receptor complexes per unit area on the cell surfaces, Pi,j(t) is the amount of ligand-releasing protease in the cell with index i,j (one index for each dimension of the two-dimensional cell lattice; for concrete indexing schemes in particular cell geometries, see Appendix), X and Y are the coordinates in the plane of the epithelial layer, Z is the transverse coordinate, and t is time. 
is the total number of ligand-receptor complexes on the surface of cell i,j. D is the diffusion coefficient of the ligand in the extracellular space, konR0 is the product of the rate of ligand-receptor forward binding and the number of receptors per unit area, koff is the ligand-receptor complex dissociation constant, and ke is the rate of ligand-induced receptor internalization.
Thus, Eq. 1 models the three-dimensional extracellular ligand diffusion in the gap of width h between the layer of cells and an impermeable barrier, whereas Eq. 2 describes the reversible ligand-receptor binding and the first-order receptor-mediated endocytosis (Lauffenburger and Linderman, 1993), Fig. 1 D. In the dynamic balance for ligand-receptor complexes, 
is the ligand concentration at the receptor-covered cell surfaces.
Similarly, Eq. 3 models the ligand-releasing protease dynamics in each cell within the epithelial layer. These dynamics consist of the combination of the first-order degradation characterized by rate constant kp and a sigmoidal generation function 
characterized by the production rate gp and the threshold CT. Such parameterization lumps a number of processes together and is common in modeling of regulatory networks (Bolouri and Davidson, 2002; Ferrell, 1997; Smolen et al., 2000). In the following we choose 
(x) to be the Heaviside function, thus assuming a sharp threshold; our results do not significantly depend on the precise form of (x). Furthermore, our first-order degradation combines the degradation at the mRNA and protein level. Let us emphasize that this thresholdlike generation term is a function of the total number of ligand-receptor complexes on the i,jth cell. Therefore, it is related to the distribution of complexes by the integral over the area of the i,jth cell surface (Eq. 3). This is supported by at least one direct measurement (Dyson and Gurdon, 1998).
The boundary condition for Eq. 1 on the surface of the epithelial layer (Z = 0), which couples the diffusion of the ligand with its secretion (as a result of intracellular processes), is given by Eq. 4. It accounts for the reversible binding and protease-mediated ligand release. The source term in this boundary condition is spatially nonuniform and varies from cell to cell across the layer. There are two contributions to this term: dissociation of ligand-receptor complexes and protease-mediated ligand release. The latter is assumed to be uniform over each cell's surface, which is expressed by the characteristic function 
i,j(X,Y), which is equal to ‘1’ on the surface of the i,jth cell and ‘0’ elsewhere, and is regulated by the availability of the ligand-releasing protease. Ligand release is modeled as first-order with respect to the protease, with the release rate per cell given by grPi,j (A is the area of the cell surface). Ligand-precursor, on which the ligand-releasing protease is acting, is assumed to be in excess. There is also a no-flux boundary condition at the impermeable barrier (Z = h) in Eq. 4.
Goals of the article and plan of the analysis
The model describes the coupled dynamics of cells in an idealized epithelial layer. To characterize these dynamics, we need to track both the extracellular and intracellular variables in the system of integrodifferential Eqs. 1–4. We would like to quantify the spatial extent of cell-to-cell communication and to analyze the effect of this coupling on the dynamics of individual cells. As an application, we consider the consequences of activating ligand release in a small group of cells within the layer. Such perturbations are implemented using techniques for tissue-specific gene expression and are routinely employed in developmental biology (Brand and Perrimon, 1993; Duffy et al., 1998). Analysis of such perturbations requires models that can resolve individual cells. It is our goal to formulate such models.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONDESCRIPTION OF THE MATHEMATICAL...RESULTSAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
characterizes the steady-state rate of ligand degradation. In this case, the spatial variation of the ligand field in the z-direction is negligible: 
This approximation is expected to be very accurate for the chosen set of geometric and dynamic parameters.
In addition to reducing the number of spatial dimensions, we reduce the number of dependent variables. Our argument is based on the timescale separation between the ligand dynamics and those of the intracellular protease. In our description, ligand binding stimulates the transcription of the ligand-releasing protease. Since this process happens on a timescale that is longer in comparison to those of binding, transport, and endocytosis, we can set the time derivatives in the dynamic balances for extracellular and receptor-bound ligand to ‘0’. This approximation imposes the following constraints on the timescales: 
and 
(see Appendix).
Together with the assumption of the ligand-limited regime, these approximations lead to the following model:
 | (5) |
 | (6) |
Hence, the problem in Eqs. 1–4 is reduced to the equation for protease dynamics in individual cells coupled to the steady linear reaction-diffusion equation for extracellular ligand. The linear problem for the ligand field can be easily solved for any particular pattern of protease activity. The resulting instantaneous field, 
can then be integrated over the area of each cell to provide arguments for the protease generation function. This leads to a fully discrete model for cell communication. In the following, this programme is implemented for square and hexagonal cells. Before that, we rescale the problem and describe the dimensionless groups.
Nondimensionalization
Eqs. 5–6 are rendered dimensionless by the following transformations,
 | (7) |
 | (8) |
|
and the maximum total number of ligand-receptor complexes per cell, 
(Tables 2 and 4).
|
|
 | (9) |
 | (10) |
and the integration is now over the rescaled cell area 
Recall that we chose 
There are only two dimensionless groups in the resulting dimensionless system
 | (11) |
To complete the derivation of the discrete model we solve the transport problem for the extracellular ligand. Since the problem is linear, the solution can be evaluated as a superposition of fields due to secretion from individual cells.
Ligand field for a single-cell source
Consider a single ligand-releasing cell placed at the origin of the epithelial layer 
, secreting the ligands with the dimensionless rate p0,0 = 1. Then, the resulting ligand field q(x,y) satisfies:
 | (12) |
everywhere on the surface of the cell placed at the origin, and ‘0’ everywhere else. The ligand fields computed for a particular set of parameters for square and hexagonal cells are shown in Fig. 2.
|
 | (13) |
,
) is the cosine transform of 0,0 (x,y). See Fig. 6, A and B, and the discussion in the Appendix. The function F(,) depends both on the shape and the size of the cell. The expressions for F(,) for the square and hexagonal cells are given in the Appendix.
|
 | (14) |
Cell-to-cell coupling coefficients
The coefficients lead to a number of useful expressions. The number of ligand-receptor complexes on the surface of the i,jth cell due to the ligand-releasing cell at the origin can be found as:
 | (15) |
is the total number of ligand-receptor complexes on the cell surface when the entire layer is producing ligand at the maximum rate. Hence, to compute the number of complexes due to a single cell, one has to multiply this maximal value 
by the interaction coefficient and divide the result by the dimensionless cell area 
: 
Using linear superposition, we express the total number of ligand receptor complexes for an arbitrary pattern of protease activity:
 | (16) |
For arrays of square and hexagonal cells these coefficients are explicitly computed in the Appendix.
Discrete model for cell communication
We now use the coupling coefficients to formulate a discrete model of cell communication. Substituting the expression for 
into the protease balance, we obtain:
 | (17) |
Parametric analysis of coupling coefficients
Given the expressions for the interaction coefficients (see Appendix), we can examine their dependence on the biophysical parameters of the original problem. In the case of EGFR signaling in Drosophila egg development, the cell area (A) and the height of the gap between the oocyte and follicle cells (h) have been estimated from microscopic images: (h 
0.5 µm, A 25 µm2; see Spradling, 1993). The extracellular diffusivity of the ligand (Spitz) and the binding/internalization rate constants can be estimated from the corresponding values in the mammalian EGFR systems. Using these parameters, we present the dependence of the coupling coefficients on the ligand diffusivity and the forward-binding rate constant in Fig. 3. Note that we plot the coupling coefficients as functions of the distance between cells. This is done only for convenience in representing discrete data points; of course in the fully discrete context, the coupling coefficients are anisotropic (although the degree of anisotropy proves to be rather small both for square and hexagonal cells). Similarly, the lines connecting the points in Fig. 3 are used only to guide the eye.
|
; Freeman, 1997; Peri et al., 2002). The rapid decay of the coupling coefficients can be exploited in the computational analysis of the discrete problem, Eq. 17. In simulating the arrays of cells, one has to evaluate the coupling only between the finite (and small) number of cells. This greatly simplifies the evaluation of the right-hand side of Eq. 17.
Analysis of the positive feedback circuit
We now illustrate the use of our model in a number of computational experiments with the positive feedback circuit, Fig. 1. Our computations are directly related to the recent results in Drosophila egg development, reviewed in Amiri and Stein (2002). In one of the experiments, Peri and co-workers used genetics to permanently activate the protease and hence, the ligand release, in a small (2–4) group of cells within the follicular epithelium (Peri et al., 2002). The ligand (Spitz) acts on the cognate receptors (EGFR) on the surfaces of epithelial cells. The authors then followed the level of expression of the gene controlled by the receptor (pipe; see Amiri and Stein, 2002). It was found that pipe was repressed both in the Spitz-releasing cells and in their neighbors. It is well known that ligand-receptor (EGFR/Spitz) binding in this system stimulates the expression of the ligand-releasing protease (Rhomboid) and that the released ligand (Spitz) can activate the protease in the neighboring cells; see Fig. 1, B and E. So what prevents the perturbation from spreading across the cellular layer? What controls the effect—the strength of ligand release, the extracellular transport, or the size of the perturbation (the number of cells constitutively expressing the protease)? Our computations illustrate how these questions might be addressed within the presented biophysical framework.
First, we construct stable stationary solutions of the discrete problem of cell-to-cell communication, Eq. 17. Let us emphasize that these solutions have no counterparts or analogs in the continuous reaction-diffusion models with bistable nonlinearity. For illustrative purposes, we consider simple three-cell perturbations in the square and hexagonal lattices, Fig. 4. In experiments, cells within these clusters expressed the protease independently of ligand-receptor binding. In terms of our model, this means that the protease level in these cells is constant. For the Heaviside nonlinearity, the stable level of the protease can have only two values: the "on" state, for which pi,j = 1, and the "off" state, where pi,j = 0.
|
on" transition (ignition) of one of the cells next to the original perturbation can be written as:
 |
In writing the "ignition" condition, we keep the dependence of coupling coefficients on 
, the dimensionless group combining the kinetic and transport properties (refer to the definition in Table 3). Using this condition, the value of the critical threshold in the right-hand side can be computed for every value of . This defines an upper boundary of the localized three-cell patterns in Fig. 4, A and B. Thus, the construction of the stability boundary amounts to evaluating only a small number of coefficients.
|
and CT (Eq. 11), we translate these dimensionless groups into the dimensional parameters of the reaction transport problem in Eqs. 1–4. In Fig. 4 we plot these boundaries as functions of ligand-release rate and the ligand-receptor affinity, the two parameters that have been shown to regulate the operation of autocrine loops in a number of cell culture EGFR systems. In our computations, the range of ligand affinity was dictated by the large amount of binding data available for the EGFR system (Lauffenburger and Linderman, 1993; Wiley et al., 2003). The rate of ligand release suggested by our computations spans the range attainable in the experiments with cultured autocrine EGFR-expressing cells (DeWitt et al., 2001; Dong et al., 1999). At this time, there are no quantitative data on ligand-release rates in vivo.
In addition to static perturbations, such as those used by Peri and co-workers (Peri et al., 2002), we can consider the perturbations that transiently activate ligand release in a group of cells. What is the outcome of a transient perturbation? Will it decay to zero or persist when the stimulus is turned off? The condition under which a pattern is extinguished (i.e., at least one of the cells within the cluster undergoes the transition to the "off" state) can also be easily formulated as a simple equation for a small group of coupling coefficients. For example,
 |
|
ibyl et al., 2003). For low rate of ligand release, the transient activation of protease activity induces only a transient response, Fig. 5 B. An interesting feature of the two-parameter diagrams in Fig. 4, A and B, is the presence of a clear minimum in the dependence of the critical rate of ligand release on the ligand-receptor binding affinity. In terms of the model, it means that the domain of parameters for which the particular localized pattern exists is disconnected. This is illustrated by a one-dimensional cut through the two-parameter diagram, Fig. 4 C. We verified that overstepping both the left and the right boundaries in this plot generates an ignition front (see the discussion above). The nature of these transitions can be described as follows. For low ligand affinities, binding and transport cannot generate the number of ligand-receptor complexes that are necessary for activating the positive feedback. On the other hand, very high binding affinities prevent an efficient transport of ligand. This regime corresponds to the case, when ligand is degraded before it is passed on to the neighbors.
Finally, Fig. 4 D demonstrates that the effect of the localized perturbation is critically affected by the size of the perturbation. For the examples presented in Fig. 4 D, the upper stability boundary that was correlated with the initiation of "ignition" fronts is negatively correlated with the size of the perturbation. The fact that the stability boundary strongly depends on the size of the perturbation is an immediate consequence of the fact that coupling beyond nearest neighbors is important. Experiments with constitutively active ligand release are frequently limited by the ability to control the size of the perturbation (Peri et al., 2002). We suggest that the corresponding results should be interpreted with care. In particular, for the same values of binding, signaling, and transport parameters, changing the size of the perturbation can move the system between the regimes of localized and long-range signaling.
CONCLUSIONS
We have developed discrete models of cell-to-cell communication in epithelial layers. In analyzing autocrine and paracrine signals, we consider the molecules, the feedback loops, and the geometry of cell-to-cell communication that are conserved across species (Casci and Freeman, 1999). EGFR system serves as a paradigm for autocrine/paracrine tissue regulation, and in multiple developmental and pathological contexts EGFR is controlled by the positive feedback loop discussed in this article (Wells, 1999). Indeed, the EGFR/Ras/MAPK-mediated feedback from ligand binding to ligand release operates in many mammalian systems (Dent et al., 1999; Doraiswamy et al., 2000; Gechtman et al., 1999; Montero et al., 2002). At this time all models of EGFR system are formulated at the level of a single cell (Wiley et al., 2003). Our work is aimed at the development of mechanistic models at the tissue level.
One of our main results is the derivation of cell-to-cell coupling coefficients as a function of geometric, cellular, and molecular parameters of the ligand transport problem. These coefficients provide a quantitative framework for the analysis of cell-to-cell interactions by diffusing ligands in epithelial layers. Our explicit expressions for the coupling coefficients obtained for particular geometries can be used to estimate the range of cell communication in epithelial layers. Note, however, that our approach—"removing" the continuum part of the problem by slaving it to intracellular variables through the introduction of cell-to-cell coupling coefficients—is not limited to periodic arrays of cells of simple shapes and can be applied to arbitrary cell arrangements.
In addition to the general analysis of autocrine signals in epithelial layers, we have analyzed a discrete and nonlinear model of the positive feedback between ligand release and binding. In particular, we found a class of stable stationary solutions in the form of clusters of a few signaling cells. The existence of these solutions is due to the essential discreteness of the considered system. These results can be applied to the EGFR signaling in Drosophila oogenesis. This is possible due to the well-characterized "geometry" of cell communication in this problem, and to the availability of molecular and cellular data for the EGFR system. The results on localized patterns can be used to analyze the patterned states generated by the localized activation of ligand release in the follicular epithelium (Pai et al., 2000; Peri et al., 2002).
While this work has been primarily motivated by the EGFR-mediated cell communication, recent advances in cell biology of developmental signaling pathways enable the formulation of mechanistic models for other systems. Recent data on ligand-receptor affinity in the Wingless pathway can be combined with the spatially resolved measurements of Wingless transport (Dubois et al., 2001; Lloyd et al., 2002; Pfeiffer et al., 2002). This provides a good incentive for the development of transport models that could account for the processes of intracellular ligand trafficking (Entchev et al., 2001; Seto et al., 2002).
The rapid decay of coupling coefficients suggested by our computations is in line with the conclusions of experimental studies of the EGFR system, both in vitro and vivo. We suggest that, in epithelial layers, many autocrine and paracrine networks may be operating in the regime of "almost" next-nearest neighbor coupling. If, as a result of future quantitative experiments, this turns out to be the case, then modeling of epithelial layers might draw from a large body of mathematical results available for lattice dynamical systems; see Cahn et al., 1998; Chow et al., 1998 for example.
Our approach has a number of limitations. While the ligand-limited regime and the "thin-fin" approximation are likely to hold for a large number of developmental contexts, special care has to be paid to assessing the validity of the assumption about the separation of the timescales corresponding to binding and transport and the intracellular variables. New methods have to be developed for the regime when this condition is not satisfied. In addition, we assumed that the generation of the ligand-releasing protease is a thresholdlike function of the current level of ligand-receptor complexes. More complex parameterizations of this dependence can be considered. We verified numerically that similar results are obtained if one uses a sufficiently sharp Hill function for the sigmoidal nonlinearity.
Another important issue that needs to be addressed is the effect of the noise on the signaling patterns due to the low number of molecules involved. In fact, a simple calculation shows that for realistic parameters (see Table 1) the number of ligand molecules over a single cell must be of order unity for the ligand-limiting regime to be realized. This must be reconciled with the use of the continuous approximation for ligand diffusion (Eq. 1). Let us point out, however, that the relevant quantity for the signaling patterns here is the number of ligand-receptor complexes per cell 
which acts as an input to the protease production. In contrast to the number of ligands, the number of complexes turns out to be large; for realistic parameters, 103. This is essentially due to the fact that the average lifetime of the ligand-receptor complexes (ke + koff)-1 is much longer than the timescale of ligand-receptor binding h/ks (see Table 1). Therefore, for these values of the effect of the fluctuations in the number of extracellular ligands will be averaged over times (ke + koff)-1 and to the leading order can be neglected. In this sense, Eq. 1 should be viewed as the equation for the probability density of finding a ligand molecule in a given infinitesimal volume element. We have performed preliminary Monte Carlo simulations that confirm these statements.
The combination of kinetic models with the detailed information about the subcellular structures has been very successful in the analysis of intracellular events (Slepchenko et al., 2002). We argue that a similar approach is critical for the analysis at the tissue level. A rapidly growing number of well characterized developmental systems makes this approach both necessary and feasible (Freeman and Gurdon, 2002).
|
APPENDIX |
|---|
|
TOPABSTRACTINTRODUCTIONDESCRIPTION OF THE MATHEMATICAL...RESULTSAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
 | (A1) |
According to this equation, for 
, the variation of s in the z-direction is negligible. Hence, 
Using this fact and averaging Eq. A1 over the vertical coordinate, we get, approximately,
 |
 | (A2) |
 | (A3) |
 | (A4) |
Eqs. A2–A4 contain four dimensionless parameters: 
For the thin-fin approximation in Eq. A1, (i.e., ratio of the geometrical and the dynamical length) has to be small. This condition is satisfied for the typical parameter set used in the model analysis (see Table 3). Furthermore, the relative timescales of extracellular ligand and ligand-receptor complexes are small: 
and 
; see Table 3. The steady-state approximation for these variables leads to Eqs. 9 and 10 (in the original scaling).
Coupling coefficients for square cells
Now we compute the coupling coefficients 
for square cells of size Lx. Let F(,), be the cosine transform of 0,0(x,y):
 | (A5) |
 | (A6) |
Coupling coefficients for hexagonal cells
For the hexagonal cells with dimensionless width 2v = 2Lv/L, the cosine transform F(,) of 0,0 is given by
 | (A7) |
is the side of the hexagon (see Fig. 6 B).
From Eq. 14, the interaction coefficient is 
where the function q(x,y) is defined by Eq. 13. To compute the coupling coefficients, we integrate this expression over the m,nth cell. This requires evaluating the integral 
With the indexing scheme (see Fig. 6 D), the position of cell centers on the hexagonal lattice is given by
 | (A8) |
for the m,nth cell, can be written as:
 | (A9) |
 | (A10) |
General comments
The integrals given by Eqs. A6 and A10 can be easily computed numerically. We have used adaptive integration routines in Mathematica and MATLAB for this purpose. The computation of a single coefficient takes several seconds on a very modest PC. Given the fact that coupling coefficients decay rather quickly as a function of the lattice indices, the vector field for the discrete problem is constructed very efficiently. In the current form, the coupling coefficients depend on the dimensionless cell size (l, or v and r) and 
a dimensionless parameter inherited from the original three-dimensional problem. By rescaling the variables of integration in the final expressions (Eqs. A6 and A10), that 
Thus, up to a constant, the integrals in A6 and A10 depend on a single dimensionless parameter 
Returning to the dimensional parameters, we find that 
is equal to the cell size, l (or v and r), normalized by the length scale for ligand variations along the surface 
The original cell communication system operates in a discrete regime when the ratio of these length scales is >1. The cells are essentially uncoupled for very large values of . Small values of correspond to the continuum regime.
|
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONDESCRIPTION OF THE MATHEMATICAL...RESULTSAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
This work was supported by National Science Foundation grant DMS-0211864.
|
FOOTNOTES |
|---|
ibyl's present address is Institute of Chemical Technology, Prague, Czech Republic. Submitted on November 26, 2002; accepted for publication February 11, 2003.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONDESCRIPTION OF THE MATHEMATICAL...RESULTSAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
Bergmann, A., M. Tugentman, B. Shilo, and H. Steller. 2002. Regulation of cell number by MAPK-dependent control of apoptosis. a mechanism for trophic survival signaling. Dev. Cell. 2:159–170.[Medline]
Bolouri, H., and E. H. Davidson. 2002. Modeling transcriptional regulatory networks. Bioessays. 24:1118–1129.[Medline]
Brand, A. H., and N. Perrimon. 1993. Targeted gene expression as a means of altering cell fates and generating dominant phenotypes. Development. 118:401–415.[Abstract]
Cahn, J. W., J. Mallet-Paret, and E. S. Van Vleck. 1998. Traveling wave solutions for systems of ODEs on a two-dimensional spatial lattice. SIAM J. Appl. Math. 59:455–493.
Carmena, A., E. Buff, M. S. Halfon, S. Gisselbrecht, F. Jimenez, M. K. Baylies, and A. M. Michelson. 2002. Reciprocal regulatory interactions between the Notch and Ras signaling pathways in the Drosophila embryonic mesoderm. Dev. Biol. 244:226–242.[Medline]
Casci, T., and M. Freeman. 1999. Control of EGF receptor signalling: lessons from fruitflies. Cancer Metast. Rev. 18:181–201.[Medline]
Chow, S. N., J. Mallet-Paret, and W. X. Shen. 1998. Traveling waves in lattice dynamical systems. J. Diff. Eq. 149:248–291.
Collier, J. R., N. A. Monk, P. K. Maini, and J. H. Lewis. 1996. Pattern formation by lateral inhibition with feedback: a mathematical model of delta-notch intercellular signalling. J. Theor. Biol. 183:429–446.[Medline]
Dent, P., D. Reardon, J. Park, G. Bowers, C. Logsdon, K. Valerie, and R. Schmidt-Ullrich. 1999. Radiation-induced release of transforming growth factor alpha activates the epidermal growth factor receptor and mitogen-activated protein kinase pathway in carcinoma cells, leading to increased proliferation and protection from radiation-induced cell death. Mol. Biol. Cell. 10:2493–2506.
DeWitt, A., J. Dong, H. Wiley, and D. Lauffenburger. 2001. Quantitative analysis of the EGF receptor autocrine system reveals cryptic regulation of cell response by ligand capture. J. Cell Sci. 114:2301–2313.[Medline]
Dong, J. Y., L. K. Opresko, P. J. Dempsey, D. A. Lauffenburger, R. J. Coffey, and H. S. Wiley. 1999. Metalloprotease-mediated ligand release regulates autocrine signaling through the epidermal growth factor receptor. Proc. Natl. Acad. Sci. USA. 96:6235–6240.
Doraiswamy, V., J. Parrot, and M. Skinner. 2000. Expression and action of transforming growth factor alpha in normal ovarian surface epithelium and ovarian cancer. Biol. Reprod. 63:789–796.
Dubois, L., M. Lecourtois, C. Alexandre, E. Hirst, and J. P. Vincent. 2001. Regulated endocytic routing modulates wingless signaling in Drosophila embryos. Cell. 105:613–624.[Medline]
Duffy, J., D. Harrison, and N. Perrimon. 1998. Identifying loci required for follicular patterning using directed mosaics. Development. 125:2263–2271.[Abstract]
Dyson, S., and J. B. Gurdon. 1998. The interpretation of position in a morphogen gradient as revealed by occupancy of activin receptors. Cell. 93:557–568.[Medline]
Entchev, E., A. Schwabedissen, and M. Gonzalez-Gaitan. 2001. Gradient formation of the TGF-beta homolog Dpp. Cell. 103:981–991.
Ferrell, J. 1997. How responses get more switch-like as you move down a protein kinase cascade. Trends Biochem. Sci. 22:288–289.[Medline]
Freeman, M. 1997. Cell determination strategies in the Drosophila eye. Development. 124:261–270.[Abstract]
Freeman, M. 2000. Feedback control of intercellular signalling in development. Nature. 408:313–319.[Medline]
Freeman, M., and J. B. Gurdon. 2002. Regulatory principles of developmental signaling. Annu. Rev. Cell Dev. Biol. 18:515–539.[Medline]
Gechtman, Z., J. L. Alonso, G. Raab, D. E. Ingber, and M. Klagsbrun. 1999. The shedding of membrane-anchored heparin-binding epidermal-like growth factor is regulated by the Raf/mitogen-activated protein kinase cascade and by cell adhesion and spreading. J. Biol. Chem. 274:28828–28835.
Hatini, V., and S. DiNardo. 2001. Divide and conquer: pattern formation in Drosophila embryonic epidermis. Trends Genet. 17:574–579.[Medline]
Hogan, B. 1999. Morph. Cell. 96:225–233.[Medline]
Hsu, T., D. McRackan, T. S. Vincent, and H. Gert de Couet. 2001. Drosophila Pin1 prolyl isomerase Dodo is a MAP kinase signal responder during oogenesis. Nat. Cell Biol. 3:538–543.[Medline]
Lander, A. D., W. Nie, and F. Y. Wan. 2002. Do morphogen gradients arise by diffusion? Dev. Cell. 2:785–796.[Medline]
Lauffenburger, D. A., and J. J. Linderman. 1993. Receptors: Models for Binding, Trafficking, and Signalling. Oxford University Press, New York.
Lloyd, T. E., R. Atkinson, M. N. Wu, Y. Zhou, G. Pennetta, and H. J. Bellen. 2002. HRS regulates endosome membrane invagination and tyrosine kinase receptor signaling in Drosophila. Cell. 108:261–269.[Medline]
Meir, E., G. von Dassow, E. Munro, and G. M. Odell. 2002. Robustness, flexibility, and the role of lateral inhibition in the neurogenic network. Curr. Biol. 12:778–786.[Medline]
Monk, N. A. M. 1998. Restricted-range gradients and travelling fronts in a model of juxtacrine cell relay. Bull. Math. Biol. 60:901–918.[Medline]
Montero, J. C., L. Yuste, E. Diaz-Rodriguez, A. Esparis-Ogando, and A. Pandiella. 2002. Mitogen-activated protein kinase-dependent and -independent routes control shedding of transmembrane growth factors through multiple secretases. Biochem. J. 363:211–221.[Medline]
Nilson, L. A., and T. Schupbach. 1999. EGF receptor signaling in Drosophila oogenesis. Curr. Top. Dev. Biol. 44:203–243.[Medline]
Othmer, H. G., and L. E. Scriven. 1971. Instability and dynamic pattern in cellular networks. J. Theor. Biol. 32:507–537.[Medline]
Owen, M. R., J. A. Sherratt, and S. A. Myers. 1999. How far can a juxtacrine signal travel? Proc. Roy. Soc. Ser. B. 266:579–585.
Owen, M. R., J. A. Sherratt, and H. J. Wearing. 2000. Lateral induction by juxtacrine signaling is a new mechanism for pattern formation. Dev. Biol. 217:54–61.[Medline]
Pai, L., G. Barcelo, and T. Schupbach. 2000. D-cbl, negative regulator of the EGFR pathway, is required for dorsoventral patterning in Drosophila oogenesis. Cell. 103:51–61.[Medline]
Peri, F., C. Bokel, and S. Roth. 1999. Local Gurken signaling and dynamic MAPK activation during Drosophila oogenesis. Mech. Dev. 81:75–88.[Medline]
Peri, F., M. Technau, and S. Roth. 2002. Mechanisms of Gurken-dependent pipe regulation and the robustness of dorsoventral patterning in Drosophila. Development. 129:2965–2975.[Medline]
Pfeiffer, S., S. Ricardo, J. B. Manneville, C. Alexandre, and J. P. Vincent. 2002. Producing cells retain and recycle Wingless in Drosophila embryos. Curr. Biol. 12:957–962.[Medline]
Pibyl, M., C. B. Muratov, and S. Y. Shvartsman. 2003. Long-range signal transmission in autocrine relays. Biophys. J. 84:883–896.
Ruohola-Baker, H., E. Grell, T. B. Chou, D. Baker, L. Jan, and Y. N. Jan. 1993. Spatially localized Rhomboid is required for establishing the dorsal-ventral axis in Drosophila oogenesis. Cell. 73:953–965.[Medline]
Sapir, A., R. Schweitzer, and B. Z. Shilo. 1998. Sequential activation of the EGF receptor pathway during Drosophila oogenesis establishes the dorsoventral axis. Development. 125:191–200.[Abstract]
Seto, E. S., H. J. Bellen, and T. E. Lloyd. 2002. When cell biology meets development: endocytic regulation of signaling pathways. Genes Dev. 16:1314–1336.
Shvartsman, S. Y., C. B. Muratov, and D. A. Lauffenburger. 2002. Modeling and computational analysis of EGF receptor-mediated cell communication in Drosophila oogenesis. Development. 129:2577–2589.
Slepchenko, B. M., J. C. Schaff, J. H. Carson, and L. M. Loew. 2002. Computational cell biology: spatiotemporal simulation of cellular events. Annu. Rev. Biophys. Biomol. Struct. 31:423–441.[Medline]
Smolen, P., D. A. Baxter, and J. H. Byrne. 2000. Modeling transcriptional control in gene networks—methods, recent results, and future directions. Bull. Math. Biol. 62:247–292.[Medline]
Spradling, A. C. 1993. Developmental Genetics of Oogenesis. The Development of Drosophila Melanogaster. Cold Spring Harbor Laboratory Press, Plainview, New York. pp. 1–70.
Urban, S., J. R. Lee, and M. Freeman. 2002. A family of Rhomboid intramembrane proteases activates all Drosophila membrane-tethered EGF ligands. EMBO J. 21:4277–4286.[Medline]
Van Buskirk, C., and T. Schupbach. 1999. Versatility in signaling: multiple responses to EGF receptor activation during Drosophila oogenesis. Trends Cell Biol. 9:1–4.[Medline]
von Dassow, G., E. Meir, E. M. Munro, and G. M. Odell. 2000. The segment polarity network is a robust developmental module. Nature. 406:188–192.[Medline]
Wasserman, J. D., and M. Freeman. 1998. An autoregulatory cascade of EGF receptor signaling patterns the Drosophila egg. Cell. 95:355–364.[Medline]
Wells, A. 1999. EGF receptor. Int. J. Biochem. Cell Biol. 31:637–643.[Medline]
Wiley, H. S., S. Y. Shvartsman, and D. A. Lauffenburger. 2003. Computational modeling of the EGF receptor system: a paradigm for systems biology. Trends Cell Biol. 13:43–50.[Medline]
This article has been cited by other articles:
 |
|
 |
|
  M. Coppey, A. M. Berezhkovskii, S. C. Sealfon, and S. Y. Shvartsman Time and Length Scales of Autocrine Signals in Three Dimensions Biophys. J., September 15, 2007; 93(6): 1917 - 1922. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 D. L. Nikolic, A. N. Boettiger, D. Bar-Sagi, J. D. Carbeck, and S. Y. Shvartsman Role of boundary conditions in an experimental model of epithelial wound healing Am J Physiol Cell Physiol, July 1, 2006; 291(1): C68 - C75. |
. (B) Solution of 
, see text for details) plotted as a function of cell-to-cell distance, k, defined as the ratio of the Euclidian distance divided by the cell size. (A) Effect of ligand-receptor affinity kon (square cells, Lx = 5 x 10-4 cm, D = 1 x 10-7 cm2 s-1). (B) Effect of ligand diffusivity D (hexagonal cells, 2Lv = 5 x 10-4 cm; kon = 0.1 nM-1 min-1). All other parameters are as in 
