Analysis of the Signal Transduction Properties of a Module of Spatial Sensing in Eukaryotic Chemotaxis

Bulletin of Mathematical Biology (2003) 65, 95–128 doi:10.1006/bulm.2002.0323

Analysis of the Signal Transduction Properties of a Module of Spatial Sensing in Eukaryotic Chemotaxis J. KRISHNAN AND P. A. IGLESIAS Department of Electrical and Computer Engineering, The Johns Hopkins University, Baltimore, MD 21218, U.S.A. The movement of cells in response to a gradient in chemical concentration—known as chemotaxis—is crucial for the proper functioning of uni- and multicellular organisms. How a cell senses the chemical concentration gradient surrounding it, and what signal is transmitted to its motion apparatus is known as gradient sensing. The ability of a cell to sense gradients persists even when the cell is immobilized (i.e., its motion apparatus is deactivated). This suggests that important features of gradient sensing can be studied in isolation, decoupling this phenomenon from the movement of the cell. A mathematical model for gradient sensing in Dictyostelium cells and neutrophils was recently proposed. This consists of an adaptation/spatial sensing module. This spatial sensing module feeds into an amplification module, magnifying the effects of the former. In this paper, we analyze the spatial sensing module in detail and examine its signal transduction properties. We examine the response of this module to several inputs of experimental and biological relevance. c 2002 Society for Mathematical Biology. Published by Elsevier Science Ltd. All

rights reserved.

1.

I NTRODUCTION

The phenomenon of chemotaxis involves the ability of a cell to sense the external concentration of chemical species and migrate towards higher concentrations of chemoattractants and away from higher concentrations of chemorepellants. Chemotaxis is crucial for the proper functioning of organisms, both unicellular and multicellular, and is believed to play important roles in morphogenesis, wound healing and tumor metastasis. There has been extensive work on chemotaxis and its underlying mechanisms in bacteria from both an experimental and a theoretical perspective, focusing more recently on the underlying signaling networks (Barkai and Leibler, 1997; Alon et al., 1999). The work on eukaryotic chemotaxis, especially on aspects related to the underlying signaling networks and signal transduction has been mainly of an experimental nature (Fisher, 1990; Parent and Devreotes, 1999). There have been few systematic modeling approaches to chemotaxis in eukaryotic systems. A conceptual framework outlined by Parent and Devreotes (1999) c 2002 Society for Mathematical Biology. Published by 0092-8240/03/010095 + 34 $35.00/0 Elsevier Science Ltd. All rights reserved.

96

J. Krishnan and P. A. Iglesias

suggests a simple mechanism for gradient sensing in such cells. The main idea is that directional sensing (which results after the ligand binds to the receptor) can be explained in terms of an interplay between two processes—one involving the rapid activation of an enzyme in proportion to local receptor occupancy and another, the slow activation of a diffusible enzyme in proportion to the average receptor occupancy. The former activates the production of PH-binding sites (representative of the output of the system) while the latter inhibits it. The interplay results in PHbinding sites being relatively localized and being present in greater proportion near the so-called leading edge. Recently, some new models have been proposed to describe various aspects of chemotaxis in such eukaryotic systems. A phenomenological model for chemotaxis in eukaryotes put forward by Meinhardt (1999) is based on an activator–inhibitor system of the Turing type which produces a stable activation pattern in response to the input (the external chemical concentration). In order for the model to be able to demonstrate sensitivity to changing inputs, the presence of a second inhibitor is postulated which dissipates the activation pattern allowing new activation patterns to form. The model of Meinhardt attempts to explain the response of a mobile cell to a chemical gradient. The biochemical interpretation of the model is centered around cytoskeletal rearrangements. It has been known experimentally that the ability of a cell to sense gradients can be decoupled from its motion. Specifically, Dictyostelium cells which have been immobilized by latrunculin (which inhibits actin polymerization) are still able to sense external gradients and provide a spatially polarized response (Parent and Devreotes, 1999). This suggests that the gradient sensing mechanism can largely be decoupled from the motion of the cell. Thus the gradientsensing aspect of chemotaxis can be studied and modeled without consideration of actin-polymerization and other aspects of cell movement and cytoskeletal rearrangement. This is important because ‘components’ of the signal transduction network responsible for gradient sensing can be studied largely in isolation—this includes controlled experiments elucidating various aspects of signal transduction responsible for gradient sensing and mathematical modeling of these networks and their underlying mechanisms. It should be noted, however, that there could, in general, be certain feedback effects of actin polymerization on gradient sensing. One important aspect of gradient sensing is adaptation: the response of the cell to a step change in homogeneous external concentration eventually returns to the original level. This feature needs to be incorporated into any gradient sensing model. A model proposed by Levchenko and Iglesias (2002) is aimed at describing gradient sensing in eukaryotic cells such as Dictyostelium and neutrophils. It is primarily aimed at describing gradient sensing of immobile cells, and any possible feedback effects on gradient sensing resulting from actin polymerization are not considered. This model incorporates the adaptation and gradient sensing properties of these cells. In addition, an amplification mechanism, believed to be necessary for cells to respond to very shallow gradients, is also incorporated.

Adaptation and Spatial Sensing

97

This model consists of two modules in series: an adaptation module, which takes as its input the external concentration, and provides a graded response to external gradients, and an amplification module which amplifies the output of the adaptation module. Both modules are comprised of biochemical networks of species, some of which could diffuse. The model is presented within a plausible biochemical framework. Another interesting model is that of Narang and coworkers (Narang et al., 2001). Their model is primarily aimed at describing a certain amplification mechanism which results in the formation of a soliton (pulse)-like static concentration profile of membrane phosphoinositides at the leading edge of the cell (location of the cell near the maximum external concentration) in response to an external gradient. This is formulated in terms of a biochemical network of species of the phosphoinositide cascade, that both react and diffuse. This postulates a rather different amplification mechanism than that presented in Levchenko and Iglesias (2002). It should be noted that this model could include feedback effects of actin polymerization on gradient sensing. The notion of splitting signal transduction networks into functional modules has received considerable attention recently (Hartwell et al., 1999). The underlying idea is that complex biological signal transduction networks can be best understood in terms of various functional modules embedded in them, and that these functional modules may be in some sense conserved from organism to organism. This paper deals with a detailed analysis of the adaptation module of the model presented in Levchenko and Iglesias (2002). This is a module which is built from ideas presented qualitatively in Parent and Devreotes (1999). In fact this module could itself be regarded as a quantification of the basic idea of Parent and Devreotes (1999), incorporating the adaptation property in addition. We analyze this module separately and in detail for various reasons. Firstly the basic properties of adaptation and gradient sensing are present in this module. Since the input to this model is the external concentration, essentially, and the output feeds into an amplification module, we can study the response of this module to various inputs of experimental and biological relevance. Further, the qualitative signal transduction properties of this module are robust to parameter variations. The structure of the adaptation module allows us to obtain explicit expressions for the response of the module to various inputs of interest. Finally, we can analyze both the given as well as other amplification modules with different mechanisms integrated with this module. Thus a thorough understanding of this module provides important insights for the overall modeling and analysis of gradient sensing, and eventually chemotaxis. The paper is organized as follows. Section 2 presents the adaptation module of the model presented in Levchenko and Iglesias (2002). It also discusses certain properties of this model. Section 3 analyzes the response of this model to various inputs: time-varying but spatially uniform, steady and inhomogeneous, and spatiotemporally varying signals. Inputs such as homogeneous oscillations

98

J. Krishnan and P. A. Iglesias

of external chemical concentration, uniform gradients and localized sources of external concentration are considered. Also considered are inputs corresponding to multiple sources, rotating localized sources, and rapidly switched external gradients. Section 4 presents the conclusions and outlines some future work.

2.

D ESCRIPTION OF THE A DAPTATION M ODULE

In this section, we present the adaptation module of the model described in Levchenko and Iglesias (2002), and examine some of its features. This module takes, as its input, the external signal S, which could be the concentration of the receptor–ligand complex, and is simply related to the external ligand concentration assuming very fast binding events relative to the time scale of the processes of interest. This signal mediates the activation of two enzymes, a locally acting ‘activator’ (of the response) enzyme, as well as an ‘inhibitor’ enzyme, which is also capable of long range diffusion. The deactivation of these enzymes is assumed to be constitutive [e.g., see Goldbeter and Koshland (1981)]. These two enzymes, in turn, mediate the activation and deactivation, respectively, of the response element. We present these facts in a little more detail. Both the activator and inhibitor enzymes are assumed to be found in active and inactive forms. Assuming that the total concentrations of these enzymes are fixed and equal to Atot and Itot , respectively, the equations governing the kinetics of the active forms of these enzymes (denoted by A and I , respectively) is given by, dA = −k−a A + ka0 S(Atot − A) dt dI = −k−i I + ki0 S(Itot − I ). dt From the assumption that A  Atot and I  Itot , these equations can be written as dA = −k−a A + ka S dt dI = −k−i I + ki S dt where ka = ka0 Atot and ki = ki0 Itot . One other species, known as the response element, also exists in two forms: the concentration of its active form is denoted by R ∗ , while that of its inactive form is denoted by R. The time variation for R ∗ is given by d R∗ = −k−r I R ∗ + kr A R dt

Adaptation and Spatial Sensing

99

where R + R ∗ = RT , a constant. This indicates the inhibition of this active form of the response element via the inhibitory enzyme, and ‘production’ (really activation) by the activating enzyme. We now consider the distributed system—a 1D periodic system of length L = 2π Ro mimicking the (thin) membrane of a circular cell of radius Ro . We nondimensionalize space by the cell radius, and time by a suitable time scale τ ; the concentrations A, I and S are nondimensionalized by reference concentrations A0 , I0 and S0 and R is nondimensionalized by RT . Once this is done, we arrive at the following set of PDEs: ∂a = −K −a a + K a s ∂t ∂i ∂ 2i = −K −i i + K i s + K d 2 ∂t ∂θ ∗ ∂r = −K −r ir ∗ + K r ar ∂t where K −a = τ k−a , K −s = τ k−s , K a = τ ka S0 /A0 , K i = τ ki S0 /I0 , K −r = τ k−r I0 , K r = τ kr A0 and r + r ∗ = 1. The dimensionless diffusion coefficient kd = k D τ/Ro2 , where k D is the diffusion coefficient of the inhibitor. Typical values of τ and Ro are τ = 5 s and Ro = 5 µm. This equation is solved in the domain θ ∈ [0, 2π ] with periodic boundary conditions. This is posed for a nonzero 2π-periodic signal input S. Hereafter we shall revert to lower case letters for the various rate constants and upper case for the variable concentrations. The concentration of the active form of the response element R ∗ is regarded as the output of this module. This module has for its starting point, the qualitative ‘model’ outlined in Parent and Devreotes (1999) which suggested that gradient detection arose from the interaction of two processes: the activation of a response by a locally acting enzyme and the deactivation of the response by an inhibitory enzyme. Here, the spatial aspect of gradient sensing entered only through the fact that the inhibitory enzyme is capable of diffusion—loosely speaking, the inhibition is based on global information regarding the external signal; the activating enzyme was ‘produced’ based essentially on local information regarding the external signal—it was assumed to be essentially nondiffusing. The model incorporates the diffusion of the inhibitory enzyme while regarding the activator as nondiffusing; the response element is also taken as nondiffusing. Note that one can consider these species to be weakly diffusing and the results are small perturbations of those presented here (the boundary conditions for these other variables are automatically satisfied without diffusion). This is the module for adaptation and gradient sensing (without the amplification) presented in Levchenko and Iglesias (2002). In addition to presenting the model skeleton, the authors also discuss the relationship between their model and the known biochemistry of eukaryotic gradient sensing based on the activation of

100

J. Krishnan and P. A. Iglesias

G-protein-associated chemokine receptors, in Dictyostelium and neutrophils. Various aspects of the biochemical signaling networks underlying gradient sensing and chemotaxis are still not well understood and still many biochemical details are missing. Regarding adaptation, it was previously assumed that adaptation occurred at the level of the receptor. This was an attractive idea, since that explained the adaptation of different processes which were G-protein activated. However recent experiments have indicated that adaptation does not occur at the level of G-protein, suggesting that adaptation occurs downstream of G-protein activation (Jin et al., 2000; Janetopoulos et al., 2001). The model was built around the fact that there is adaptation of the phosphoinositide phosphates P I (3, 4)P2 and P I (3, 4, 5)P3 . Based on known biochemical facts, it was suggested that the activator in the above model is phosphatidylinositol 3-kinase-γ (PI3K), while R and R ∗ were P I (4, 5)P2 and P I (3, 4, 5)P3 . Also a likely candidate for the inhibitor was suggested to be PTEN. The experiments on Dictyostelium were performed to analyze its response to an external gradient of cAMP. In the model, a tentative downstream amplification mechanism, amplifying the output of the adaptation module was suggested. Various details related to the biochemistry of eukaryotic gradient sensing around which the model is built, are described in Levchenko and Iglesias (2002). We briefly mention a few features regarding the amplification mechanism presented in Levchenko and Iglesias (2002). It is recognized that the output of the adaptation module is further amplified, and a tentative amplification mechanism was suggested. The amplification module employed involves a positive feedback mechanism. Here the reaction product R ∗ is involved in supplying a reactant for a new reaction, though it is not involved in the activation of the reaction. This new reaction is part of a network involving the positive feedback. If the external gradient is removed, supply of the substrate for this reaction is stopped and the positive feedback is deactivated. An important feature of this module is the fact that a continuous sensitivity to the external input is maintained, as a result of the nature of the amplification. The amplification module exhibits a switch-like behavior, which implies that input signals above a certain threshold are amplified considerably, while those below the threshold are not amplified. 2.1. Adaptation property. We shall analyze the above module in some detail. One important property which this model possesses is that of adaptation. This means that when the input S is homogeneous, the eventual output (in this case R ∗ ) returns to the basal level, following a step change in S. That this is the case, is seen easily as follows. Since the solution has no spatial or temporal dependence, we have, I = (ki /k−i )S A = (ka /k−a )S R ∗ = kr A/(k−r I + kr A).

Adaptation and Spatial Sensing

101

This last expression simplifies to R ∗ = (A/I )/(k−r /kr + A/I ). Since A/I does not depend on the signal S (when it is homogeneous), it immediately follows that neither does R ∗ . Note that this feature does not depend on particular values of the kinetic rate constants and is thus robust to parameter variation. The value of R ∗ corresponding to spatially homogeneous inputs is given by ∗ Radap =

ka k−i /k−a ki . k−r /kr + ka k−i /k−a ki

(1)

This is a property of the kinetics of the enzymes. Adaptation is a special feature of signal transduction of this network. From the perspective of control engineering, this suggests that a constant (step change) input does not affect the eventual output. This is known in the control engineering literature as a step rejection response. It is known that if a nonlinear system is able to robustly reject step inputs, then its linearization about a particular state will exhibit the presence of an integral control mechanism (Csete and Doyle, 2002). Details about this fact are discussed in Francis and Wonham (1975) and Byrnes and Isidori (1990). The consequences for bacterial chemotaxis are discussed in Yi et al. (2000). We present a few more details on this fact. A common way to stabilize a dynamic system at a particular state (which may be unstable) is to apply some form of feedback control. In the simplest case, known as proportional feedback, the control involves applying an effort which is proportional to the deviation of the current state of the system from the desired state (called the set point). For a simple one degree of freedom system, this results in a modified dynamic system du = f (u) + K (u − u 0 ) dt where the second term is the proportional feedback employed to stabilize the variable u at the value u 0 . K is called the gain of the feedback, and is in general a matrix. It is chosen so as to stabilize the state u = u 0 . Note that u − u 0 is the error between the state and the desired state. Another control mechanism involves applying an effort which is proportional to the integral of the error. This results in a dynamic system of the form du = f (u) + K dt

t

Z 0

(u(t) − u 0 ) dt

102

J. Krishnan and P. A. Iglesias

which can be written as du = f (u) + K y dt dy = u(t) − u 0 . dt Similarly a proportional derivative control involves applying an effort which is proportional to the time derivative of the deviation of a variable from its desired value. Often in control engineering practice, a combination of proportional and integral control (PI control) or a combination of proportional, integral and derivative control (PID control) is employed. We consider the kinetics of our adaptation module, and linearize about a particular state I0 , A0 , R0 corresponding ˆ Rˆ and S, ˆ the to an external signal S0 . Denoting the deviation variables by Iˆ, A, linearized equations are written as d Iˆ = −k−i Iˆ + ki Sˆ dt d Aˆ = −k−a Aˆ + ka Sˆ dt d Rˆ ˆ 0 + kr A. ˆ = −(k−r I0 + kr A0 ) Rˆ − (k−r Iˆ + kr A)R dt The presence of an integral control mechanism can be seen by considering the ˆ where evolution of a variable k1 Iˆ + k2 Aˆ + k3 R, k1 =

k−r R0 k−i (k−r I0 + kr A0 )

k2 =

−ki k−r R0 ka k−i (k−r I0 + kr A0 )

k3 =

−1 . k−r I0 + kr A0

Denoting this variable by Zˆ , we have d Zˆ ˆ = R. dt This demonstrates the explicit presence of an integral control action which is seen ˆ Iˆ, Zˆ variables. Thus more clearly when the full system is written in terms of A, ˆ any steady state of the full system has to satisfy R = 0, which is just the adaptation property. It should also be noted that in our case this steady state is guaranteed to be stable.

Adaptation and Spatial Sensing

103

We also briefly discuss the inclusion of the nonlinear terms. For a given state characterized by the quantities S0 , I0 , A0 , R0∗ , if we consider the variable Zˆ as described before, and include the nonlinear terms, we obtain the following equation for the time rate of change of Zˆ :   d Zˆ k−r Iˆ + kr Aˆ ˆ = R 1+ . dt k−r I0 + kr A0 This can be written in a simplified form as   d Zˆ k I + k A −r r = Rˆ . dt k−r I0 + kr A0 Since the quantity in brackets is always nonzero, for any nonzero signal, it follows that at steady state Rˆ equals zero. It may also be mentioned that the existence of an integral feedback action can also be demonstrated in other direct ways. It should be noted that arguments of the kind presented above are based solely on specific signal transduction properties of a network and robustness. The argument asserts that if a certain kind of special property (e.g., rejection of a particular class of input signals) is possessed by the network and if this feature is robust, then the mathematical structure of the dynamics of the network has a special feature; further this can be interpreted in terms of a control mechanism embedded in it. Note, however, that because of the nature of this argument it is not often possible to attribute simple biological interpretations to these control mechanisms, particularly in complex networks. Nevertheless, it may be possible by a procedure similar to that above, to be able to identify a subset of species in a network responsible for such special signal transduction properties.

3.

R ESULTS

3.1. Homogeneous inputs. We consider the response of the module under consideration, to various kinds of external inputs. Most of the inputs considered here have some experimental or biological relevance. As mentioned earlier, immobilizing the Dictyostelium cells using latrunculin makes it possible to study aspects of gradient sensing independent of motion. This also makes it possible to study aspects of the underlying networks experimentally in a more systematic and careful manner, and interpret the results more easily than in the case of mobile cells. Such experiments are being currently performed by different groups. In these experiments immobilized cells are subjected to different inputs and the response of the cell is observed. As a measure of the response of the cell, the fluorescence along the membrane is tracked: green fluorescent protein tagged CRAC molecules move from the cytosol and attach to certain PH binding sites on the membrane.

104

J. Krishnan and P. A. Iglesias

While this is related to the membrane concentration of PIP3(3,4,5), the correlation is not quite so clear since there are other phosphoinositides to which the CRAC molecules bind (Servant et al., 2000). Thus we are currently mainly concerned with qualitative features of the signal transduction. We note that immobilizing the cell prevents any feedback effects of actin polymerization on the gradient sensing network, so that this additional feature would have to be taken into account when modeling gradient sensing in mobile cells. The simplest kind of external input is a uniform steady concentration and this leads to a uniform steady output as discussed above. A step increase in the (homogeneous) external concentration causes the concentration of both the activator and the inhibitor to increase, while that of the (active form of) response element first increases, reaches a maximum and then decreases (as a result of the inhibition effect) back to its original value (correspondingly similar behavior is seen when a step decrease is applied). The time taken for the response element to reach its (transient) maximum can be easily determined in terms of the various rate constants. For all other plots in this paper the following parameters are used: ki = 1.0, k−i = 1.0, ka = 4.0, k−a = 2.0, kr = 3.0, k−r = 3.0 and kd = 0.70. This dimensionless diffusion coefficient corresponds to a diffusion coefficient of 3.5 µm2 s−1 . Since the values of various kinetic constants are not known, we have chosen a particular set for illustration purposes. These are chosen taking into account a few features: the activation and deactivation of the activator enzyme A is faster than the corresponding processes for the inhibitor; also the rate constants for the activation and deactivation of the response element are more or less the same. Of course, one of the critical issues in trying to build and analyze models of biological networks is to determine how robust they are to parameter variations. We will perform a detailed analysis of the quantitative dependence of the response on parameter values. Analytical solutions provide a transparent description of the response in terms of inputs and system parameters. We consider the case of a homogeneous oscillatory input. In the simplest case, we apply a ‘square-wave’ type input consisting of periodic step increases in external concentration followed by step decreases of equal height and duration. If the time period of this square-wave is much greater than the time scale of the various reactions in our module, then the response will be one of transient increase of response element concentration followed by a recovery back to the original state during the first part of every cycle, and a transient decrease of response element concentration followed by a recovery back to the original state in the second part of every cycle, resulting in a periodic response. This is seen in Fig. 1 where periodic steps of height h = 0.5 (relative to the average level) are applied as input; the time duration for each step is T = 8. As is seen in the figure, a steady oscillatory response of the kind described above is obtained. Note also that the deviations of response element concentration from the steady value are not the same in the positive and negative parts of the cycle.

Adaptation and Spatial Sensing

105

Figure 1. Evolution of activator (dashed line) and inhibitor (solid line) concentrations (left), and response element concentration (right) in response to a rectangular wave of step height 0.5 and total period 16.

We consider analytically the response of this module to various oscillatory inputs: periodic inputs where each period comprises of a positive step followed by a negative step of equal height and of equal duration, and periodic sinusoidal inputs, all superimposed on some steady uniform external concentration. The response in all cases is uniform and oscillatory. The response is easily determined by solving the equations with an oscillatory input. The input S = a + b f (t), where f (t) is periodic with period T . The inhibitor equation reads dI + k−i I = ki S(t). dt The solution is easily obtained as I (t) = e

−k−i t

t

Z

ki ek−i τ S(τ ) dτ + I (0) e−k−i t .

0

This shows that for an oscillatory input S(t), the inhibitor concentration is oscillatory and the dependence on the initial condition I (0) decays with time. The oscillatory solution which is approached is proportional to ki . The effect of increasing k−i is to decrease the inhibitor concentration in general: increasing k−i causes the inhibitor concentration at any fixed time instant to be reduced. The solution for the activator equation is obtained just as above, and again the same conclusions with regard to the solution and its dependence on the forward and backward rate constants hold. Z t −k−a t A(t) = e ka ek−a τ S(τ ) dτ + A(0) e−k−a t . 0

The solution of the equation for the response element is also obtained in a similar way. The equation for the concentration of active form of the response element is

106

J. Krishnan and P. A. Iglesias

Figure 2. Evolution of activator (dashed line) and inhibitor (solid line) concentrations (left), and response element concentration (right) in response to a sinusoidal input S(t) = 0.6 + 0.5 sin(2πt/16).

written as

d R∗ + (k−r I (t) + kr A(t))R ∗ = kr A(t) dt

where I (t) and A(t) are obtained from the previous equations. Again, the solution is similar to that above: R ∗ (t) = e−

Rt

0 (k−r I (u)+kr

A(u)) du

t

Z

kr A(τ ) e

Rτ 0

(k−r I (y)+kr A(y)) dy

dτ

0

+ R ∗ (0)e−

Rt

0 (k−r I (u)+kr

A(u)) du

.

Again, we see that the dependence on the initial condition decays with time, and stable oscillations are observed. Here too, increasing the deactivation constant k−r causes an overall decrease in R ∗ . Increasing the production constant kr increases the production of R ∗ . The dependence of the amplitude of oscillations on the various problem parameters can also be determined. The response of the module to a sinusoidal input of the same total period T = 16 and amplitude h = 0.5 as before was considered. The results are shown in Fig. 2. The system again asymptotes to steady oscillations. However in the sinusoidal case, the total amplitude of the response element concentration is clearly less than in the rectangular wave cases. This is because the system is not exposed for enough time to the maximum concentration of the signal S as in the rectangular wave case. It may be noted that a steady oscillatory response was also seen in response to a rectangular wave input consisting of steps of equal height but unequal duration. Another point worth mentioning in this context is that while the eventual response of the module adapts when subject to a step response, that is not the case when the input is a periodic function like a sine wave. For the robust rejection of a sine wave, the module would have to contain in some form a model of the

Adaptation and Spatial Sensing

107

signal (with the same frequency) to be rejected, along with feedback (Francis and Wonham, 1975). It should also be noted that there exist nongeneric parameter values for which the system can reject oscillatory inputs. 3.2. Spatial sensing. The previous subsection dealt with the response of the model to various homogeneous but time varying inputs. The response depended only on the kinetic features of the model. However the main interest in this module lies in its spatial sensing properties. In this subsection, we consider spatially varying but static inputs. The response in these cases will be spatially varying but steady, and the response will depend on the diffusion properties of the inhibitor, in addition to the kinetic properties of the network. We start by considering simple inputs like a constant gradient, and then consider other more complex inputs. We start by considering a situation corresponding to an external concentration gradient of constant strength, in the x direction. The concentration of the chemical species outside the cell, can be written as S = a 0 + b0 x (assuming negligible distortion from the cell) where x is the dimensionless x-coordinate, which in turn can be written as S = a 0 + b0 Ro cos θ on the surface of the cell, shaped as a disk of (dimensionless) radius Ro centered around the origin. Thus the dependence may be written as S(θ ) = a + b cos θ . We now determine the steady-state response of our module to this input. The steady-state distribution is obtained simply as A(θ) =

ka (a + b cos θ ) k−a

and is proportional to the input. The concentration profile of the inhibitor is obtained by solving the steady-state equations: kd

d2 I − k−i I + ki (a + b cos θ ) = 0 dθ 2

with periodic boundary conditions. This is easily solved using a Fourier series decomposition of the inhibitor profile. The inhibitor profile has only the zeroth and first cosine modes and is given by I (θ ) =

ki ki a+ b cos θ. k−i k−i + kd

Thus the inhibitor profile displays a similar variation as the input, a cosine variation, but has less relative amplitude because of diffusion. The steady-state profile of the active form of the response element is obtained simply as R ∗ (θ ) =

A(θ )/I (θ ) . k−r /kr + A(θ )/I (θ )

108

J. Krishnan and P. A. Iglesias

Note that the factor A(θ )/I (θ ) =

(ka /k−a )(a + b cos θ) (ki /k−i )(a + b cos θ/(1 + α))

where α = kd /k−i . We note several features of this response. Both the activator and inhibitor have cosine variations in the angle, just like the external input; thus the maxima and minima of the activator and inhibitor coincide with that of the input. By analyzing the expression of the response element profile, we see that the extrema of the response element coincides with that of the activator and inhibitor concentration profile. Further, the maximum of the activator/inhibitor profile, also corresponds to a maximum of the response element profile (for any positive diffusion coefficient kd ). Thus the response element concentration profile, in this simple case, has the features which we expect: its maximum is located at the location of the maximum of the external signal on the cell surface. We note that the value of R ∗ depends only on the relative gradient b/a of the external concentration. Thus the imposition of a constant external gradient results in a persistent (steady and stable) spatially varying response. This response of course depends on the various kinetic parameters; in addition it also depends on the diffusivity of the inhibitor. One characteristic of the response which is usually measured experimentally is the ratio of the response concentration at the front of the cell (i.e., where the concentration is maximum) to that at the back of the cell (where the concentration is minimum). This ratio is determined from above as R ∗ (front) =

(K a /K i )(a + b)/(a + b/(1 + α)) K r + (K a /K i )(a + b)/(a + b/(1 + α))

R ∗ (back) =

(K a /K i )(a − b)/(a − b/(1 + α)) K r + (K a /K i )(a − b)/(a − b/(1 + α))

R ∗ (front)/R ∗ (back) =

(K a /K i )(a + b)/(a + b/(1 + α))(K r + (K a /K i )(a − b)/(a − b/(1 + α)) (K a /K i )(a − b)/(a − b/(1 + α))(K r + (K a /K i )(a + b)/(a + b/(1 + α)))

where in the above equation, K i = ki /k−i , K a = ka /k−a , K r = k−r /kr . Thus the dependence of this quantity on various parameters is explicitly known and easily computed. Figure 3 shows the activator, inhibitor and response element concentrations, when subject to an input S = 1.0 + 0.4 cos θ clearly demonstrating the features just discussed. Note that the variation of the gradient will result in a response which is qualitatively as shown, but with a different amplitude. We consider the effect of the various parameters on the system response. While the qualitative features of the response do not change, quantitative features naturally change. In the above equations we see that response element characteristics depend on the external signal through the ratio b/a which is the relative gradient.

Adaptation and Spatial Sensing

109

Figure 3. Steady-state activator (dashed line) and inhibitor (solid line) concentration profiles (left), and response element concentration profile (right) in response to an input S(θ) = 1.0 + 0.4 cos θ , corresponding to a linear gradient.

Also, while the response depends on various kinetic parameters, the dependence in basically through two groups of parameters: K = K a /(K i K r ) and α. Note that the group K is just a combination of kinetic parameters and occurred even when one considered the adaptation phenomenon—the basal level response depends on this group of parameters. The other group α is where the diffusion coefficient of the inhibitor is present. A detailed study of the effect of variation of parameters is performed in the Appendix. The response element concentration profile, when subjected to other steady, but inhomogeneous inputs is also easily obtained just as above: A(θ ) = (ka /k−a )S(θ ) R ∗ (θ ) =

A(θ )/I (θ ) . k−r /kr + A(θ )/I (θ)

The inhibitor profile is obtained by expanding it as a Fourier series and doing the same for the external input: I (θ) =

∞ X (ck cos kθ + dk sin kθ) k=0

S(θ ) =

∞ X ( pk cos kθ + qk sin kθ). k=0

The coefficients ck , dk are obtained from the inhibitor equation directly as ck = ki pk /(k−i + kd k 2 ) dk = ki qk /(k−i + kd k 2 ).

110

J. Krishnan and P. A. Iglesias

Figure 4. Input profile corresponding to one localized source (left), and two symmetrically placed localized sources (right), shown as a function of angle.

Thus the inhibitor profile, and the response element profile can be reconstructed for any input. Besides the constant gradient case, there are other inputs of interest. One such input is a localized ‘pulse-like’ input, which mimics the presence of a localized source of concentration, such as a needle used experimentally (Parent and Devreotes, 1999). The concentration profile of a single pulse-like input and that corresponding to two symmetric pulse-like inputs are shown in Fig. 4. The response to these inputs was also obtained numerically using a standard finite difference scheme in MATLAB. Figure 5 shows the steady-state behavior of the module when subjected to the single pulse input shown previously. We again see that the response element concentration in response to a localized input has a maximum near the location of the maximum of the external concentration profile. We note that experimentally a relatively localized response near the maximum of the external concentration is seen. Existing models such as Narang et al. (2001) and Meinhardt (1999) do predict a localized structure near the maximum of the external concentration; in the case of Meinhardt’s model, a clear dependence on the external signal cannot be made since his model is really an amplification mechanism. Narang’s model predicts a soliton-like structure which does not depend on the features of the external signal (gradient and mean value), as long as the external gradient is strong enough to form the soliton. One of the interesting issues in chemotaxis, is the response of the system to multiple chemotactic cues. This could, in general, be due to multiple sources of the same chemical, or be a result of the response to more than one chemical. A detailed understanding of the chemotactic response to multiple cues is very important for understanding chemotaxis in biologically realistic situations. Another reason for studying the response of the cell to multiple sources is that it sheds light on the internal interaction and integration of various responses to external stimuli. This is especially true in the context of gradient sensing in immobilized cells where these issues can be studied experimentally in a systematic and controlled manner.

Adaptation and Spatial Sensing

111

Figure 5. Steady-state activator (dashed line) and inhibitor (solid line) concentration profiles (left), and response element concentration profile (right) in response to the ‘single pulse’ input shown previously.

Some experiments on the response of chemotaxing cells to more than one source of concentration of the same chemical have been performed (Foxman et al., 1997, 1999). This was done for mobile cells, and the results were not very conclusive. Even the simpler case of the response of immobile cells to multiple sources of concentration of the same chemical has not been systematically examined experimentally. A systematic study of the response of immobile cells to multiple sources of the same chemical, would provide valuable insights into the signal processing capabilities of the chemotactic network. We perform a couple of calculations where there are multiple sources of external concentration. First we consider the case of perfectly symmetrically placed localized inputs, at diametrically opposite points. Thus the input concentration S(θ ) has only even Fourier modes, and in addition has a profile which has two spikes. This profile was shown in Fig. 4(b). The solution of the equations governing the activator, inhibitor and response elements results in a profile for the response element (and also the activator and inhibitor) which also possesses this symmetry. In addition, this profile is stable. Thus the response of the adaptation module to such an input is to provide a response which is equally poised towards initiating motion towards each of the two sources. The activator, inhibitor and response element concentration profiles for this case are shown in Fig. 6. We note that this response is independent of how the two sources are ‘assembled’ or their history. We also note that a stable symmetric response is always obtained irrespective of parameter values. However, if the sources are not of equal strength, the response of the module will depend on the relative strengths of the sources. In this case, the response element concentration profile will have two local maxima, localized near each of the sources, but a slightly higher maximum near the stronger of the sources. Such an unequal double-pulse input is shown in Fig. 7(a), where the pulses are symmetrically located, but of different amplitude. Figure 8 shows the steady-state

112

J. Krishnan and P. A. Iglesias

Figure 6. Steady-state activator (dashed line) and inhibitor (solid line) concentration profiles (left), and response element concentration profile (right) in response to the symmetric double pulse input shown previously. An ‘equal’ and symmetric response is obtained.

Figure 7. Input profiles corresponding to an unequal double pulse (left) and a double humped pulse (right).

activator, inhibitor and response element profiles, showing clearly some bias for the ‘stronger’ source. It is clearly seen, that now the response is in favor of the stronger source, but exhibits local maxima near both locations of the needles. We also note that this does not depend on how the two sources were introduced (i.e., simultaneously or consecutively). The general conclusion is that this part of the gradient sensing network will not exclusively bias the response towards any one source. We note that the model of Narang will cause a symmetric response even in this case, since the eventual response does not depend on the detailed nature of the sources as long as the input is strong enough to excite the system (locally) past a threshold. Thus his model will predict two identical solitons at the two local maxima corresponding to the (unequal) sources. Yet another possibility is the case where there are two localized sources placed close to each other. The response of an immobilized cell to such an input, could

Adaptation and Spatial Sensing

113

Figure 8. Activator and inhibitor concentration profiles, and response element concentration profile in response to the unequal double pulse input shown previously.

Figure 9. Steady-state activator (dashed line) and inhibitor (solid line) concentration profiles (left), and response element concentration profile (right) in response to an input corresponding to two localized sources relatively close to each other shown previously.

indicate certain important facts of the gradient sensing network, especially with regard to the amplification mechanism. An experimental scenario related to this would involve placing the two ‘needles’ (localized sources of external chemical concentration) close to each other. Such an input is shown in 7(b). The response of the module corresponding to an input of two pulses relatively close to each other is shown in Fig. 9(b). We see that the profile of the inhibitor is rather ‘smoothed out’ as a result of diffusion; however the profile of the response element exhibits two maxima indicating that it is still able to distinguish between these sources. It would be of much interest to determine experimentally how the system would respond in this case, as this will suggest to what extent the response is able to distinguish between two closely placed sources. Overall, the results in this subsection indicate that the response to steady inhomogeneous inputs is a steady and persistent inhomogeneous output. The presence

114

J. Krishnan and P. A. Iglesias

of multiple sources of the same chemical results in a response which ‘sees’ the multiple inputs but is biased towards the strongest source. 3.3. Spatiotemporal inputs. In this section we consider the response of the adaptation module to inputs with both spatial and temporal dependence. We have previously seen that the module is capable of responding to purely temporal variations of external concentration, as well as static inhomogeneous inputs. It is also important to examine if the gradient sensing property persists when the location of an external source is changed and if it is able to show the same sensitivity as before. We consider three cases where the location of an external source is changed—in two of the cases, the external source is changed continuously, and in one case, it is changed abruptly. The first input is a rotating pulse—a localized peak of external concentration which rotates with an angular velocity ω. This mimics a rotating needle, an input which has been used experimentally (Devreotes, 2001). Since the input has time dependence, so will the concentration profiles of the activator, inhibitor and response element. However, since the rotation is steady we expect, and will see that this is indeed the case, that the response is also steadily rotating in time with the same angular velocity. This means that in a frame rotating with angular velocity ω, the profiles of activator, inhibitor and response element appear steady. Just as before, the original governing equations are ∂I ∂2 I = −k−i I + ki P(θ − ωt) + kd 2 ∂t ∂θ ∂A = −k−a A + ka P(θ − ωt) ∂t ∂ R∗ = −(k−r I + kr A)R ∗ + kr A ∂t where P(θ ) is a pulse-like concentration profile, similar to that seen in the previous section [Fig. 4(a)]. By looking for steadily rotating solutions, I (θ, t) = i(θ − ωt), A(θ, t) = a(θ − ωt) and R ∗ (θ, t) = r ∗ (θ − ωt), we find that the steadily rotating profiles satisfy the following ODEs (with independent variable φ = θ − ωt): kd

∂ 2i ∂i +ω − k−i i + ki P(φ) = 0 2 ∂φ ∂φ ω ω

∂a − k−a a + ka P(φ) = 0 ∂φ

∂r ∗ − (k−r i + kr a)r ∗ + kr a = 0. ∂φ

We can now solve this steady problem, just as before. Expanding the input P and

Adaptation and Spatial Sensing

115

the activator and inhibitor profiles in Fourier series, we have, P(φ) =

∞ X

pk cos kφ + qk sin kφ

k=0

i(φ) =

∞ X

ck cos kφ + dk sin kφ

k=0

a(φ) =

∞ X

ek cos kφ + f k sin kφ.

k=0

Substituting this into the governing equation results in the equation for the Fourier coefficients of i and a. ck (kd k 2 + k−i ) − ωkdk = ki pk dk (kd k 2 + k−i ) + ωkck = ki qk . Thus the Fourier coefficients of the inhibitor profile are given by ck = ki ( pk (kd k 2 + k−i ) + ωkqk )/((kd k 2 + k−i )2 + ω2 k 2 ) dk = ki (− pk ωk + qk (kd k 2 + k−i ))/((kd k 2 + k−i )2 + ω2 k 2 ). Those of the activator profile are 2 + ω2 k 2 ) ek = ka ( pk k−a + ωkqk )/(k−a 2 f k = ka (− pk ωk + qk k−a )/(k−a + ω2 k 2 ).

From this, the entire activator and inhibitor profile can be reconstructed. The response element profile can be obtained by directly solving the corresponding differential equation dr ∗ − (1/ω)(k−r i(φ) + kr a(φ))r ∗ = −kr a(φ)/ω dφ where a(φ) and i(φ) are known. Denoting −(1/ω)(k−r i(φ) + kr a(φ)) by F(φ) and −kr a(φ)/ω by G(φ), the solution to this equation is r ∗ (φ) e

Rφ 0

F(ψ) dψ

φ

Z

G(ψ) e

=

Rψ 0

F(δ) dδ

dψ + r ∗ (0)

0

r ∗ (φ) = e−

Rφ 0

F(ψ) dψ

φ

Z

G(ψ) e 0

Rψ 0

F(δ)dδ

dψ + r ∗ (0) e−

Rφ 0

F(ψ) dψ

.

116

J. Krishnan and P. A. Iglesias

Figure 10. Activator (dashed line) and inhibitor (solid line) concentration profiles (left), and response element concentration profile (right) in the rotating frame in response to a rotating pulse input with ω = 0.1. Also shown in dotted lines are these profiles for the case when ω = 0.5.

This last equation can be solved for the unknown constant r ∗ (0) by requiring that r ∗ (0) = r ∗ (2π ). This gives the complete response element profile. Thus the entire profile (a rotating wave: I = i(θ −ωt), A = a(θ −ωt), R ∗ = r (θ −ωt)) is obtained using the procedure above. Figure 10 shows the profiles of activator, inhibitor and response element in the rotating frame in response to a rotating ‘pulse-like’ profile of the same shape as shown before (Fig. 4). The result can be compared with the corresponding responses to the stationary input of the same shape. It may also be mentioned that experiments performed, involving a rotating needle do indeed show a rotating relatively localized response (Devreotes, 2001) (also see some images of preliminary experiments at www. dictybase.org/tutorial/rotatingwave). It should be emphasized that these experiments will involve the effects of amplification mechanisms too, which are not discussed here. Figure 10(a) and (b) show the concentration profiles of the activator, inhibitor and response element (in the rotating frame), in response to this rotating pulse input, for two values of the rotation speed ω = 0.1 and 0.5. What one observes here is that for the higher ω value, the peak of the response lags the peak of the input profile (which is at θ/2π = 0.5) to a greater extent. It would be of interest to compare the response of models with other amplification mechanisms for this experiment. We now turn to the analysis of the module response to an input similar to what has been briefly described in Postma and Van Haastert (2001): fairly rapidly switching input gradients. Thus a constant input gradient is imposed in one direction (say the x-direction) and after a steady response has been established, the direction of the gradient is quickly reversed keeping the mean concentration fixed. This is done until the new steady response is established, and a periodic switching of gradients is imposed in this manner. One way to effect this kind of gradient change is also to rotate the (immobile) cell by 180◦ periodically in a fixed concentration

Adaptation and Spatial Sensing

117

field (in this case, the spatial coordinates used above are those in a frame attached to the cell). We can consider both the rapid as well as a smoother oscillatory transition. Thus we model the input S (for an immobile cell placed at the origin) as S = a + b cos θ H (t) where H (t) is either a smooth oscillation, e.g., H (t) = sin ωt or a more abrupt transition such as a rectangular wave (switching between H = 1 and H = −1), of total period T . We thus consider the equations ∂I ∂2 I = −k−i I + ki (a + b cos θ H (t)) + kd 2 ∂t ∂θ ∂A = −k−a A + ka (a + b cos θ H (t)) ∂t ∂ R∗ = −(k−r I + kr A)R ∗ + kr A. ∂t We solve these equations in a straightforward manner similar to that before. We first decompose the activator and inhibitor profiles into Fourier spatial modes: I = I0 + I1 cos θ A = A0 + A1 cos θ. Note that I0 , I1 , A0 , A1 are all functions of time. The governing equations for the quantities I0 , I1 , A0 , A1 are d I0 dt d I1 dt d A0 dt d A1 dt

= −k−i I0 + ki a = −(k−i + kd )I1 + ki bH (t) = −k−a A0 + ka a = −k−a A1 + ka bH (t)

whose solution is I0 (t) = e−k−i t

t

Z

ki ek−i τ adτ + I0 (0) e−k−i t

0

I1 (t) = e

−(k−i −kd )t

t

Z

ki e(k−i +kd )τ H (τ ) dτ + I1 (0) e−(k−i −kd )t

0

A0 (t) = e−k−a t

t

Z

ka ek−a τ adτ + A0 (0) e−k−a t

0

A1 (t) = e−k−a t

t

Z 0

ka ek−a τ H (τ ) dτ + A1 (0) e−k−a t .

118

J. Krishnan and P. A. Iglesias

We see thus, that aysmptotically I0 → ki a/k−i , a constant, while I1 (t) asymptotically approaches an oscillatory solution, with period T , just as seen before. Exactly the same conclusions follow for the functions A0 (t) and A1 (t). This determines the profiles of the inhibitor and activator I (t) and A(t), and shows that they asymptotically approach a (inhomogeneous) state which oscillates in time, with a fixed mean value. They oscillate between a state where the maximum is at one end θ = 0 and one where the maximum is at the other end θ = π. Finally the response element concentration can be obtained by solving the corresponding equation: Z t Rt Rτ ∗ − 0 (k−r I (θ,u)+kr A(θ,u)) du R (θ, t) = e kr A(θ, τ ) e 0 (k−r I (θ,y)+kr A(θ,y))dy dτ 0

+ R ∗ (θ, 0) e−

Rt

0 (k−r I (θ,u)+kr

A(θ,u)) du

.

This solution reveals the simple fact that the response element concentration profile is also oscillatory, and varies from one state where the maximum of the profile is at one end (θ = 0) to another state where the maximum of the profile is at the other end (θ = π). Figure 11 shows the space–time plot of the response element concentration for this case. In the space–time plot darker regions correspond to higher concentration and lighter regions to lower concentrations. We see various ‘bands’ where the maximum of the profile is either at θ = 0 (the edges of the plot) or at θ = π. We also see a rapid transition between these adjacent bands. We clearly see that the maximum of the profile changes from end to end, with every switching. This is also demonstrated in Fig. 12 which shows the time variation of the response element concentrations at the leading (θ = 0) and trailing (θ = π ) edges. Thus the maximum of the response element concentration profile is able to ‘relocate’ to the location of the new maximum of the input concentration profile. We also perform simulations for the case where a localized source similar to that shown earlier is kept at one end (θ = π/2), and rapidly removed and placed at an opposite end (θ = 3π/2). Equivalently the immobile cell is rotated through 180◦ . Qualitatively exactly the same scenario described above is seen in this case also. This aspect is shown in Fig. 13(a). A small point worth noting is that in addition to the response element having a maximum which develops at the location of the new maximum, it has a smaller ‘bump’ at the location of the old maximum. This in turn is due to the inhibitor profile having a small depression at this location. The response element concentration profile after having adjusted to the external concentration profile in a new cycle is shown in Fig. 13(b). The results of the experiments of Van Haastert have not been described in any detail. However, the ‘relocation’ of the response to the position of the new maximum every cycle is mentioned and this apparently happens in a robust manner. We finally consider the response of the adaptation module to a traveling wave input. A traveling wave input is of much biological importance, since in the process of aggregation, Dictyostelium cells signal to each other in the form of traveling waves. There is a fairly extensive literature on the signaling wave patterns that are

Adaptation and Spatial Sensing

119

Figure 11. Response element concentration profile as a function of space and time in response to the periodically switched linear gradient. This plot clearly shows the location of the maxima changing—for one part of the cycle, the maximum is at θ = 0 and for the other part of the cycle, the maximum is at θ = π. Higher concentrations are denoted by pink regions and lower concentrations by violet regions in this plot. The time period of switching is T = 8.

1 0.9 0.8 0.7 s

0.6 0.5 0.4 0.3 0.2 0.1 0 0

2

4

6

8

10 12 14 16 18 20 x

Figure 14. Traveling wave concentration profile S(z) (one wavelength is shown); space– time plot of concentration of response element in response to the traveling wave input— darker regions indicate higher concentrations.

responsible for intercellular communication. For recent references dealing with chemotaxis in such a context and the extracellular signaling patterns which form, see Goldstein (1996), Palsson and Cox (1996), Palsson et al. (1997) and Lauzeral et al. (1997). We consider an input which is a traveling wave in the x-direction. We make the simplifying assumption that the cell does not significantly distort this traveling wave pattern. Strictly speaking, the traveling wave has to go ‘around’ the cell creating a genuine two-dimensional pattern. Given that the equations governing the extracellular wave propagation are not easily solved analytically

120

J. Krishnan and P. A. Iglesias

Figure 12. Response element concentration at θ = 0 (solid line) and θ = π (dashed line) as a function of time. It is clearly seen how the location of the maximum of the concentration profile, for one part of the cycle, becomes the location of the minimum of the concentration profile for the other.

(a)

(b)

Figure 13. Response of the cell to rapidly switching the location of an external localized source similar to that shown earlier. (a) Time variation of variation of response element concentration at the ‘front’ (θ = π/2, solid line) and ‘back’ (θ = 3π/2, dashed line) of the cell (the locations at which the external concentration is alternatingly maximum). Note that the maximum ‘relocates’ to the location of the new external maximum. (b) Response element concentration profile, showing a small ‘bump’ at the location of the previous maximum. This is also seen in (a), where the concentration of the response element on the side (dotted line) is lower than that of both the back and the front.

we make this simplification. The main purpose of performing this calculation, with the additional approximation, is to get some qualitative insight. A more detailed calculation of the extracellular field in specific cases can also be performed (Palsson and Cox, 1996; Dallon and Othmer, 1998). The traveling wave input is denoted by W (x − ct), where c is the speed of the wave; this is a function which is periodic with period L. This implies that the input S at any location is time-periodic with period L/c. Thus the chemotactic

Adaptation and Spatial Sensing

121

Figure 15. Activator and inhibitor concentrations (left), and response element concentration (right) at the locations θ = 0 (solid line) and θ = π (dashed line) in response to the traveling wave input shown previously. Note that the inhibitor and response element concentrations oscillate with different amplitudes at the two locations.

input to the cell is periodic in θ (with period 2π) and also periodic in time with period L/c. Thus the eventual response of the cell—the concentrations of the activator, inhibitor and response element also have exactly the same feature—is periodic in space and time. The asymptotic response of the cell can be easily derived in terms of the input. The derivation of the relevant expressions is described in the Appendix. The response of the module to a traveling wave of profile shown in Fig. 14(a) is considered. This traveling wave has speed c = −2 traveling in the −x direction and wavelength L = 20 which corresponds to 10 cell diameters. The response element concentration as a function of space and time is shown in Fig. 14(b). We see that at any spatial location, the response element concentration is oscillatory in time. But looking at Fig. 15(a) and (b), we see that the amplitudes of the oscillations vary. In particular we see that the amplitude of the response element concentration is higher at θ = 0, the ‘edge’ which sees the wave first, when compared to θ = π. Looking at Fig. 15(b) we see that while the activator concentration oscillates with the same amplitude at both locations, the inhibitor concentration oscillates with a higher amplitude at θ = π, due to the building up of inhibitor here (due to diffusion) even before the wave ‘reaches’ this location. As a result the amplitude of oscillation of the response element concentration at the ‘back’ of the cell is less than that at the ‘front’. Of course a realistic solution of this problem would have to consider the effect of the cell too, and this feature is likely to further decrease the amplitude of the response element oscillation at the location θ = π . 4.

C ONCLUSIONS

In this paper, we have analyzed in detail the signal transduction properties of the adaptation module of the model proposed in Levchenko and Iglesias (2002)

122

J. Krishnan and P. A. Iglesias

for spatial sensing of eukaryotic cells such as Dictyostelium and neutrophils. The model presented in Levchenko and Iglesias (2002) is aimed at describing the main features of gradient sensing in these cells—adaptation to homogeneous external signals, an inhomogeneous and persistent response to inhomogeneous steady signals and the ability to sense shallow external concentration gradients. The model consists of two modules: a so-called adaptation module which takes as its input a signal which is simply related to the external concentration (for example, the concentration of a receptor–ligand complex, which is related to the external concentration under a quasi-steady state hypothesis), and an amplification module which takes as its input, the output of the adaptation module. The adaptation module itself can be regarded as a quantification and extension of the conceptual framework of gradient sensing outlined in Parent and Devreotes (1999). The output of the adaptation module is the concentration of the active form of a species referred to as a response element. In order to test and validate the model, a detailed analysis of this model and its signal transduction properties is necessary. This paper dealt with a detailed analysis of the adaptation module. A detailed analysis of the amplification module presented in Levchenko and Iglesias (2002) and comparison with other amplification mechanisms will be considered in a future publication. The adaptation module, while being fairly simple, is still very important since it possesses the basic adaptation and spatial sensing properties. The response of this module to various kinds of inputs—temporal and homogeneous, steady and inhomogeneous, and spatiotemporal is necessary to understand the working of this module, and generates predictions which can be tested. Most of the inputs considered had experimental/biological relevance. The sensitivity of the response to problem parameters is very important, as it indicates qualitative trends of variation of the response to the change of the rates of various reactions in the underlying network. Such trends could be compared to those obtained experimentally, by promoting or inhibiting any of the reactions by genetic or biochemical means. It should be mentioned that the qualitative features of the signal transduction properties of this network as presented here are robust to parameter variation. While the model demonstrated adaptation to homogeneous steady signals, the response to homogeneous oscillatory signals is oscillatory. The module also produces an inhomogeneous steady and persistent response to steady external inhomogeneous signals. For a constant gradient, the response exhibits a maximum exactly at the location where the external concentration was maximum. Other cases, like a localized spike in external concentration, were also studied. The response of the module to multiple sources of the chemical input was also examined. This indicates that the module is capable of ‘seeing’ and distinguishing between different inputs, and is in general biased towards the strongest source. In the case of spatiotemporal input signals, stable spatiotemporal responses were observed. These results indicate that the gradient-sensing properties are valid for time-varying external gradients. All these studies were performed using analytical

Adaptation and Spatial Sensing

123

and numerical methods. In many of the cases considered, explicit expressions for the response of the module to inputs was obtained. Presently we are analyzing the amplification module in Levchenko and Iglesias (2002) as well as other biochemically based amplification mechanisms. Since nonlinearity appears to be an important component of any amplification mechanism, the nature of the nonlinearity could strongly affect the qualitative features of signal transduction. A systematic analysis of various properties of these amplification mechanisms could suggest qualitative differences in signal transduction, which could in turn suggest experiments, which would discriminate between them. This would give us a better picture of gradient sensing, which eventually could be integrated with cell movement to result in a more detailed understanding of chemotaxis in eukaryotic cells.

ACKNOWLEDGEMENTS We thank Peter Devreotes and Andre Levchenko for useful discussions. We gratefully acknowledge financial support from the Whitaker foundation and NSF under grant DMS-0083500. A PPENDIX Derivation of traveling wave response. The response of the module under consideration to traveling wave inputs is presented here. Suppose the traveling wave is denoted by W (x − ct) where c is the speed of the wave. W is a periodic function of x − ct with period L, the wavelength of the wave. Since the same holds good for the function S (with x now being the x-coordinate of a point on the cell surface), we can expand S in a Fourier series. Denoting z = x − ct, we first decompose S into Fourier modes:

S(z) =

∞ X

G m cos 2πmz/L + Hm sin 2π mz/L .

m=0

The various constants G m and Hm depend on the traveling wave profile. Since z = Ro cosθ − ct (in our case Ro = 1 based on the nondimensionalization performed), this expression can now be written as

S(θ, t) =

∞ X m=0

∞ X

 Gm

2πm Ro cos θ 2π mct 2πm Ro cos θ 2πmct cos cos + sin sin L L L L



  2π m Ro cos θ 2πmct 2π m Ro cos θ 2πmct + Hm sin cos − cos sin . L L L L m=0

124

J. Krishnan and P. A. Iglesias

This can be written, in simplified notation as S(θ, t) =

∞ X

gm (θ) cos

m=0

2π mct 2πmct + h m (θ) sin L L

where the functions gm and h m are 2π periodic functions of θ and are easily obtained from above. We can now further expand the functions gm and h m in Fourier series in θ . This results in an expansion of the form S(θ, t) =

∞ X ∞ X

em,n cos

m=0 n=0

+

∞ X ∞ X

2πmct 2π mct cos nθ + f m,n cos sin nθ L L

gm,n sin

m=0 n=0

2π mct 2πmct cos nθ + h m,n sin sin nθ L L

which is just the double Fourier expansion of S in space and time: it is periodic in θ with period 2π and periodic in time with period L/c, the time period of the traveling wave. To obtain the (asymptotic) activator and inhibitor profiles, we expand these profiles in similar double Fourier series, anticipating that the asymptotic state is one which is periodic in time with period L/c. I (θ, t) =

∞ X ∞ X

am,n cos

m=0 n=0

+

∞ X ∞ X

2πmct 2π mct cos nθ + bm,n cos sin nθ L L

cm,n sin

m=0 n=0

A(θ, t) =

∞ X ∞ X

pm,n cos

m=0 n=0

+

∞ X ∞ X m=0 n=0

2πmct 2πmct cos nθ + dm,n sin sin nθ L L

2πmct 2πmct cos nθ + qm,n cos sin nθ L L

rm,n sin

2πmct 2πmct cos nθ + sm,n sin sin nθ. L L

Plugging in to the governing equations, we obtain algebraic equations for the various coefficients. For the inhibitor profile, these are m(2πc/L)cm,n + (k−i + kd n 2 )am,n = ki em,n −m(2πc/L)am,n + (k−i + kd n 2 )cm,n = ki gm,n m(2πc/L)dm,n + (k−i + kd n 2 )bm,n = ki f m,n −m(2πc/L)bm,n + (k−i + kd n 2 )dm,n = ki h m,n .

Adaptation and Spatial Sensing

125

The first two equations are linear equations (for each m and n) from which one can solve for am,n and cm,n in terms of em,n and gm,n . Similarly, the next two equations can be solved for bm,n and dm,n in terms of known quantities f m,n and h m,n . This allows one to construct the asymptotic, oscillatory spatiotemporal profile of the inhibitor concentration. The case of the activator follows in exactly the same way. The equations are m(2πc/L)rm,n + (k−a ) pm,n = ka em,n −m(2πc/L) pm,n + (k−a )rm,n = ka gm,n m(2πc/L)sm,n + (k−a )qm,n = ka f m,n −m(2πc/L)qm,n + (k−a )sm,n = ka h m,n . Thus the asymptotic spatiotemporal profile of the activator is obtained. Both the activator and inhibitor profiles are periodic in time (with period L/c) in addition to being periodic in space. To obtain the asymptotic profile of the response element, we solve the response element equation exactly as before. Since the asymptotic profiles of I and A are periodic in time, we solve the equation for R ∗ as an ODE in time, for every θ, with coefficients which are periodic in time. This provides the asymptotic solution for R ∗ (which can also be obtained numerically). Z t Rt Rτ ∗ − 0 (k−r I (θ,u)+kr A(θ,u)) du R (θ, t) = e kr A(θ, τ )e 0 (k−r I (θ,y)+kr A(θ,y)) dy dτ 0

+ R ∗ (θ, 0)e−

Rt

0 (k−r I (θ,u)+kr

A(θ,u)) du

.

Effects of parameter variation. While the qualitative behavior of the module is robust to parameter variation, there is a quantitative effect and that is what we discuss briefly here. We will study the effect of parameter variation in the context of a linear external concentration, which results in S(θ) = a + b cos θ . As indicated in the text, the response is persistent, with maxima and minima at the location of the external maxima and minima, respectively. First we address the question of the effect of the external signal. We will use the ratio of the response element concentrations at the front and back of the cell as a measure of the gradient sensing capacity of the cell. It is easily seen from the expression for the response element concentration profile (see text) that this ratio (and for that matter even the response element concentration at any location) depends only on the relative gradient b/a of the external profile, and the variation is monotonic as expected. The more interesting case is the effect of the model parameters. The expression for the quantity of interest, as derived in the text can be rewritten as R ∗ (front)/R ∗ (back) =

K (a + b)/(a + b/(1 + α))(1 + K (a − b)/(a − b/(1 + α)) K (a − b)/(a − b/(1 + α))(1 + K (a + b)/(a + b/(1 + α)))

126

J. Krishnan and P. A. Iglesias (a)

(b)

Figure 16. Dependence of gradient sensing capabilities, as measured by ratio of response element concentrations at the front and back of the cell on (a) parameter α which is a dimensionless measure of the diffusion effects of the inhibitor (for the case K = 2) and (b) ratio of rate constants K (for the case α = 0.7).

where in these equations, K = K a /(K i K r ). It is thus seen that dependence of this quantity (and more generally the response element concentration profile) on the various parameters is in terms of two parameter groups: K and α = kd /k−i . The parameter values used for study in the text corresponded to K = 2 and α = 0.7. We first consider the effect of the variation of the inhibitor diffusion coefficient, which is described by varying α. Note that when α = 0, the case of a nondiffusing inhibitor, the ratio of the response element concentration values (hereby referred to as R f /Rb) is just unity. Thus because of the structure of the kinetics of the activator and inhibitor variation (which is what guaranteed adaptation), a nondiffusing inhibitor implies that the response element concentration is uniform everywhere. This outlines the need for inhibitor diffusion, in this context, to achieve gradient sensing. We see, in Fig. 16(a) that there is a monotonic variation of R f /Rb on the parameter α. This means greater inhibitor diffusion results in a greater spatial variation of the response element profile. We also note that in the limit of large α, the quantity R f /Rb approaches a limiting value. This limiting value is given by R f /Rb(α → ∞) =

(a + b)(1 + K (1 − b/a)) . (a − b)(1 + K (1 + b/a))

The limiting case occurring for large inhibitor diffusion, corresponds to a case where the inhibitor variation is smoothed out due to diffusion, and the response element spatial variation is basically the result of a position dependent activation, via the activator. Next, we consider the effect of variation of the parameter K —this is depicted in Fig. 16(b). We notice that the variation is again monotonic, with the gradient sensing effect decreasing as a function of increasing K . We notice that as K becomes very large, the ratio R f /Rb approaches unity. Thus, making activator ‘production’ more efficient, or inhibitor ‘production’ weaker, will decrease the

Adaptation and Spatial Sensing

127

sharpness of gradient sensing. At the other end, we note that when the constant K approaches zero, the ratio R f /Rb approaches a constant value given by R f /Rb(K → 0) =

(a + b)(a − b/(1 + α)) . (a − b)(a + b/(1 + α))

R EFERENCES Alon, U., M. G. Surette, N. Barkai and S. Leibler (1999). Robustness in bacterial chemotaxis. Nature 397, 168–171. Barkai, N. and S. Leibler (1997). Robustness in simple biochemical networks. Nature 387, 913–917. Byrnes, C. I. and A. Isidori (1990). Output regulation of nonlinear systems. IEEE Trans. Autom. Control 35, 131–140. Csete, M. E. and J. C. Doyle (2002). Reverse engineering of biological complexity. Science 295, 1664–1669. Dallon, J. C. and H. G. Othmer (1998). A continuum analysis of the chemotactic signal seen by Dictyostelium discoideum. J. Theor. Biol. 194, 461–483. Devreotes, P. N. (2001). Personal communication. Fisher, P. R. (1990). Pseudopodium activation and inhibition signals in chemotaxis by Dictyostelium discoideum amoebae. Semin. Cell Biol. 1, 87–97. Foxman, E. F., J. J. Campbell and E. C. Butcher (1997). Multistep navigation and the combinatorial control of leukocyte chemotaxis. J. Cell Biol. 139, 1349–1360. Foxman, E. F., E. J. Kunkel and E. C. Butcher (1999). Integrating conflicting chemotactic signals: the role of memory in leukocyte navigation. J. Cell Biol. 147, 577–877. Francis, B. A. and W. M. Wonham (1975). The internal model principle of control theory. Automatica 12, 457–465. Goldbeter, A. and D. E. Koshland Jr (1981). An amplified sensitivity arising from covalent modification in biological systems. Proc. Natl. Acad. Sci. 78, 6840–6844. Goldstein, R. E. (1996). Traveling wave chemotaxis. Phys. Rev. Lett. 77, 775–778. Hartwell, L. H., J. J. Hopfield, S. Leibler and A. W. Murray (1999). From molecular to modular cell biology. Nature 402, C47–C52. Janetopoulos, C., T. Jin and P. N. Devreotes (2001). Receptor mediated activation of heterotrimeric G-proteins in living cells. Science 291, 2408–2411. Jin, T., N. Zhang, Y. Long, C. A. Parent and P. N. Devreotes (2000). Dynamic localization of the G-protein βγ complex in living cells during chemotaxis. Science 287, 1034–1036. Lauzeral, J., J. Halloy and A. Goldbeter (1997). Desynchronization of cells on the developmental path triggers formation of spiral waves of cAMP during Dictyostelium aggregation. Proc. Natl. Acad. Sci. 94, 9153–9158. Levchenko, A. and P. A. Iglesias (2002). Models of eukaryotic gradient sensing: application to chemotaxis of amoebae and neutrophils. Biophys. J. 82, 50–63. Meinhardt, H. (1999). Orientation of chemotactic cells and growth cones: models and mechanisms. J. Cell Sci. 112, 2867–2874.

128

J. Krishnan and P. A. Iglesias

Narang, A., K. K. Subramanian and D. A Lauffenburger (2001). A mathematical model for chemoattractant gradient sensing based on receptor-regulated membrane phospholipid signaling dynamics. Ann. Biomed. Eng. 29, 677–691. Palsson, E. and E. C. Cox (1996). Origin and evolution of circular waves and spirals in Dictyostelium discoideum territories. Proc. Natl. Acad. Sci. 93, 1151–1155. Palsson, E., K. J. Lee, R. E. Goldstein, J. Franke, R. H. Kessin and E. C. Cox (1997). Selection for spiral waves in the social amoebae Dictyostelium. Proc. Natl. Acad. Sci. 94, 13719–13723. Parent, C. A. and P. N. Devreotes (1999). A cell’s sense of direction. Science 284, 765–770. Postma, M. and P. J. M. Van Haastert (2001). A diffusion translocation model for gradient sensing by eukaryotic cells. Biophys. J. 81, 1314–1323. Servant, G., O. D. Weiner, P. Herzmark, T. Balla, J. W. Sedat and H. R. Bourne (2000). Polarization of chemoattractor receptor signaling during neutrophil chemotaxis. Science 287, 1037–1040. Yi, T.-M., Y. Huang, M. I. Simon and J. C. Doyle (2000). Robust perfect adaptation in bacterial chemotaxis through integral feedback control. Proc. Natl. Acad. Sci. 97, 4649–4653. Received 23 April 2002 and accepted 30 September 2002