Diffusion Models for Chemotaxis: A Statistical Analysis of Noninteractive
Unicellular Movement
JOSEPH C. WATKINS Department of Mathematics, University of Southern California, Los Angeles, California 90089-1113 AND BIRGIT WOESSNER Fachbereich Mathematik, Universitiit Urn, Federal Republic of Germany Received 7 September 1989; revised 30 November 1990
ABSTRACT A program is developed for applying stochastic differential equations to models for chemotaxis. First a few of the experimental and theoretical models for chemotaxis both for swimming bacteria and for cells migrating along a substrate are reviewed. In physical and biological models of deterministic systems, finite difference equations are often replaced by a limiting differential equation in order to take advantage of the ease in the use of calculus. A similar but more intricate methodology is developed here for stochastic models for chemotaxis. This exposition is possible because recent work in probability theory gives ease in the use of the stochastic calculus for diffusions and broad applicability in the convergence of stochastic difference equations to a stochastic differential equation. Stochastic differential equations suggest useful data for the model and provide statistical tests. We begin with phenomenological considerations as we analyze a onedimensional model proposed by Boyarsky, Noble, and Peterson in their study of human granulocytes. In this context, a theoretical model consists in identifying which diffusion best approximates a model for cell movement based upon theoretical considerations of cell phsyiology. Such a diffusion approximation theorem is presented along with discussion of the relationship between autocovariance and persistence. Both the stochastic calculus and the diffusion approximation theorem are described in one dimension. Finally, these tools are extended to multidimensional models and applied to a threedimensional experimental setup of spherical symmetry.
1.
INTRODUCTION
In 1905, Einstein [ll] gave his explanation for Brownian motion via a kinetic theory of molecules in part to convince the vitalists once and for all of the hopelessness of their point of view. In a strange twist of irony scarcely a decade after Einstein’s description of the subject, Przibaum [27] MATHEMATICAL
BIOSCIENCES
OElsevier Science Publishing 655 Avenue of the Americas,
104:271-303
(1991)
Co., Inc., 1991 New York, NY 10010
271 0025-5564/91/$03.50
272
JOSEPH C. WATKINS
AND BIRGIT
WOESSNER
found that Brownian movement was an adequate description for the movement of unicellular organisms. (See Thompson [31], page 76). Seventy years later, biologists have much more technologically advanced tools than a microscope, a pencil, and a metronome to investigate the motion of such organisms. At the same time, mathematicians have much more technologically advanced tools to understand Brownian motion and related random processes. The purpose of this paper is to inject some of these probabilistic tools into the biological literature in such a way that biologists can avail themselves of them. In conformity with recent mathematical models, we consider two types of cellular movement: the swimming motion of bacteria and the migrating motion of cells moving along a substrate. These models can be used to try to understand cell movement from the more superficial phenomenological point of view. In other words, the experimenter can take data and construct a consistent model with no real regard for the underlying cellular or molecular processes that might give rise to the activities under study. More fundamentally, we can propose that the randomness in the mathematical model derives from intrinsic fluctuations in the local environment and in the sensing mechanism. Cell movement is thus modeled as a stochastic process whose randomness derives from these two sources. In the particular case of Escherichiu coli, a detailed statistical description of the movement was presented by Berg and Brown [3]. By building a microscope that automatically follows the individual cells, they found an alternating sequence of two types of behavior quite different in character. In the first the change in direction is gradual, and in the second the change in direction is abrupt. They called these two behavior types runs and twiddles. During a run the flagellar filaments of a bacterium work together in a bundle as it moves. When the cell begins to twiddle, the bundle loosens. The bundle then re-forms, and the cell travels in a new direction that depends upon the change in orientation of the bundle relative to the body of the cell. Berg and Brown observed that the change in direction between the run preceding and the one following a twiddle is an angle whose distribution does not vary from twiddle to twiddle. Also, the speed of the bacterium is essentially constant during a run and between two different runs. Twiddle times are short compared to run times, and in an isotropic environment both the twiddle times and the run times have an exponential distribution. If an attractant is introduced into the environment, the distribution of twiddle times remains unchanged. The run times still have an exponential distribution. However, the parameter in the distribution depends on the angle between the run and the gradient of the concentration of the attractant. The bacterium seems to have no control during a twiddle. Consequently, if a cell has a preferred direction, it can influence its movement only by altering the length of the run times.
DIFFUSION
MODELS
FOR CHEMOTAXIS
273
Polymorphonuclear leukocyte or granulocyte cells will serve as an example of a cell whose typical behavior is locomotion along a surface. In actuality, these cells thrust themselves forward by means of a lamellipodium. The forming of each lamellipodium involves the streaming of cytoplasm, which brings the bulk of the cell forward. Trinkhaus [33] describes this mode of locomotion in great detail and speculates on the mechanisms that might explain its behavior. A comprehensive statistical description of leukocyte movement will be quite difficult. Under optimal conditions, Zigmond et al. [36] describe the leukocyte as morphologically and behaviorally polarized. Movement is smooth and is governed by a stable lamellipodium. Under suboptimal or artificial conditions, leukocytes can make erratic changes in polarity due to competition among multiple lamellipodia (see Gerisch and Keller 1131 and Shields and Haston [291X How does one turn these observations into a mathematical model for chemotaxis? As usual, one wants to simplify the situation to produce a model that captures the essential features of the actual phenomena in a way that can be treated mathematically. Ideally, the model should be susceptible to statistical tests that check its reliability. To this date, most models of swimming bacteria assume that a run is linear motion at a constant speed and a twiddle is a change in direction. In addition, models typically assume that all cells have the same pattern of behavior and that the cells do not interact, for example via contact or chemical transmitters. When experimenters try to conduct experiments in an environment that is as homogeneous (in time or in space) as possible, the model reflects this attempt by assuming a suitably homogeneous environment. More detailed models include the exponential distribution for the random lengths of twiddle times and run times. The usual models following from these assumptions are random walks, or more generally Markov processes. The Berg-Brown experiment is prototypical for unicellular movement. Data are collected for short-time behavior and are synthesized to supply the parameters in the model. The statistical checks for the reliability of the model involve its ability to predict long-term behavior. Typical questions might include: Do the cells have a preferred direction in the presence or the absence of an attractant? What is the probability that the cell finds the attractant? How much time, on the average, does the cell require to find the attractant? Peterson and Noble [26] use a two-dimensional random walk model to analyze human granulocyte movement. If this random walk has no preferred direction-interpreted to mean that the direction of the run has mean zero-then the variance in the position of a cell should be linearly proportional to the number of runs. Peterson and Noble observed the movement of granulocytes in the absence of bacteria and subsequently in the presence of bacteria. Their data were the number of runs and the
274
JOSEPH C. WATKINS
AND BIRGIT
WOESSNER
displacement during each of the runs. The total distance traveled was computed, and the results were analyzed by performing a t-test to determine whether the mean direction of a run was indeed zero. In the absence of bacteria, the difference between the actual distance and the predicted distance was insignificant (P > O.OS>,whereas when bacteria were present, the difference was highly significant (P < 0.001). Their conclusion: The runs have mean displacement zero if no bacteria are present. However, if bacteria are present the data are extremely unlikely to be a sample of random displacements having mean zero. They interpret this unlikeliness as evidence for chemotaxis. Nossal and Weiss [21] propose a descriptive theory for cell movement in two dimensions. Assuming that the starting point for the cell is the origin of coordinates and that the attractant is located on the x axis, they derive expressions for the mean and variance of the cell displacement after n runs. Their model has considerable detail: Run and twiddle times are exponential in the manner described by Berg and Brown, the turning angle is uniformly distributed, and path lengths and turning angles are independent. Using the words of Berg and Brown, the model has an abundance of Poisson statistics, and consequently the model is essentially Markovian. A random process is called Markovian if the probability distribution of any future occurrence depends on the history of the process up to the present time only via the present situation. The exponential distribution is a typical feature in Markov processes because the probability per unit time of the termination of a twiddle or the termination of a run is constant. The past, in this case the length of the present run or twiddle, is no help in predicting when the run or twiddle will end. Thus the pair (X(t),V(t)) taken together form a Markov process. If we see the cell at X(t) moving with velocity V(t), we know that it will continue in this direction for an exponential length of time and that the parameter for this exponential random variable depends on X(t) and V(t). As we begin to make predictions about the future after time t, we see that we need only X(t) and V(t) to make these predictions. To date, Wolfgang Alt [l] has developed the most comprehensive axiomatic model for chemotaxis. His model begins with a general setup for the two-parameter stochastic process (X(t), V(t)) as a time inhomogeneous semi-Markov process. As we have seen in a Markov process, V(t) is constant for an exponential length of time. The parameter in this exponential distribution can depend on position and on time. A semi-Markov process maintains the Markov property for the sequence of velocities v,, vz, v,, . . . that a cell happens to choose, but allows the run time to take on an arbitrary distribution. As before, these run times are independent and have an identical distribution, which can depend on the location of the cell and on time. (See Othmer et al. [24] for more on this. They call this
DIFFUSION MODELS FOR CHEMOTAXIS
275
process a kangaroo process or a position jump process. See Cinlar [81 for an extensive development of semi-Markov processes.) Alt soon restricts his attention to time-homogeneousprocesses. This acknowledges the attempts that experimenters make to create a time-homogeneous environment and the difficulties in statistical inference for time-inhomogeneous processes but ignores the fact that cells fatigue or age. From biochemical considerations, Ah can further restrict his model to give plausible expressions for the sequence of velocities and for the times for which the cell maintains these velocities. For more, see Okubo [23]. In this paper, we adopt a point of view that is part of the modern trend in probabilistic analysis; that is, we focus on the random processes themselves rather than on the partial differential operators, difference equations, or integral kernels that are featured in more analytical treatments. The hope is that by taking this point of view we can find a more unified framework for the present discussion and a more natural path for generalization in future discussions. After all, the objects of study are random processes and thus they rightfully belong at the center of attention. Let us emphasize that this treatment supplements, but cannot replace, the sound analytical treatment that, for example, Wolfgang Alt has given. Indeed, Alt’s analysis is essential in identifying parameters in the model and in estimating the errors involved in using a diffusion approximation. The model advocated in this manuscript is a diffusion model. The questions that a biologist should ponder before adopting a diffusion model are two:’ (1) Is the motion continuous? (2) Does the distribution of future motion of the cell depend only on its present position? Typically, we can easily decide on an answer to the first question. An affirmative answer to question 2 is always rendered with some reservations. For example, we have just argued that direction is also needed to predict future movement in the Nossal-Weiss model. However, one may reason that bacteria make many runs, with a direction change at the end of each run. So the current direction, which the cell follows for only a very short time, has very little impact on the cell’s ability to find its attractant. If this rationale is convincing, then, for the purpose of modeling, yes is a reasonable response to question 2. Certainly, in the model described above, this response is mathematically well-grounded. We are asking a question concerning the result of a large number of small random fluctuations, and in such an instance the limit theorems of probability theory play a key role. Answering yes for question 2 is similar to accepting the sample mean of random variables as the mean and then ascribing a normal distribution about the mean. In other words, a
276 yes
JOSEPH C. WATKINS AND BIRGIT WOESSNER
to question 2 reflects the fact that the law of large numbers and the central limit theorem can be used to give good approximations. In the mathematical literature, the Nossal-Weiss model is called to transport process (see Chandresekar [7]), and answering the two questions in the affirmative would constitute an affirmation of the applicability of a version of the central limit theorem appropriate in this context. Because the details of the velocity term vanish in the limit and are replaced by average values for the velocity, the limit theorem is sometimes called homogenization. In a homogeneous environment this average quantity appears as the slope of a linearly increasing variance. In following the tracks of paramecia, Przibaum witnessed this behavior “with great exactitude.” On the other hand, this linearity was not found with great exactitude by Gruler and Biiltman 1141in their study of leukocytes. However, linear growth in the variance of displacement occurred after a period of sublinear growth. This is one of many aspects of a behavior properly called persistence. As we shall see, if the ratio between the typical time of events of interest and this persistence time is sufficiently large, then a diffusion model is appropriate. Consequently, if the variance is not asymptotically linear, then a diffusion model is never appropriate. For example, if the variance grows at some power of time, we may choose to investigate a fractal model. The definition of a diffusion is a very simple statement: A diffusion is a Markov process whose motion is continuous. Not all diffusions can be easily estimated statistically. However, if we can closely approximate the motion from each point by a Brownian motion that may not drift in some direction, then, as we shall learn, the diffusion process may be represented as the solution to a stochastic differential equation. Chemotaxis, the physical process of interest, is a multidimensional process. The theory we shall describe extends easily to higher dimensions and is complicated only by the additional notation necessary in multidimensional calculus. Although the diffusion equation appears frequently in the literature on chemotaxis, there has been, to date, no systematic exploitation of the probabilistic perspective. Here, we formulate the problem in terms of a stochastic differential equation. This formulation is neither more nor less than the classical diffusion equation but serves to introduce us to a new set of mathematical tools. For example, Boyarsky et al. [5] mention the use of stochastic differential equations. However, in the absence of a probabilistic description of the events, they cannot devise good statistical tests to evaluate their data on boundary-crossing probabilities. These boundarycrossing probabilities provide only one of several possible statistical and experimental methods upon which one may decide on a diffusion model. For example, the mean time until this boundary crossing is easily computed and provides more evidence for accepting or rejecting any given model. If the experiment can be run in several geometries, then we can be more
DIFFUSION
MODELS
FOR CHEMOTAXIS
277
certain that our data are revealing something intrinsic concerning cell movement rather than some artifact of the geometry. Thus in Section 4 we give an example with spherical symmetry to demonstrate the adaptability of the method. With this case study, it is hoped that the experimenter can establish the appropriate stochastic differential equation and the corresponding statistical tests. 2.
ANALYSIS
OF A ONE-DIMENSIONAL
MODEL
We apply these ideas to the situation in which a stochastic differential equation model was first proposed. In [S], Boyarsky et al. observed the motion of human granulocytes. Ninety to 120 cells were placed in a capillary tube, and their movement was recorded by photographing a microscopic image (54 X) every 8 s. The x component of the displacement was recorded at 80 s, the shortest time in which movement could be reliably measured. Each photograph was divided into 2.5 strips 1 cm wide, and all cells were identified with their original strip. This experiment was reproduced six times in a neutral case (control) and six times with bacteria present on the left side of the capillary tube (chemotaxis). These data are summarized in Figure 1 with histograms. A move of less than one-half the cell diameter (3 mm under magnification) is difficult to detect and hence is set to zero on the histogram. A second set of data recorded the probability of reaching zero, the left side, before 2.5, the right side of the photograph at later times. The data recorded after an experiment had run 96 min in chemotaxis and 80 min in the control case are included in Tables 1 and 2, respectively. Histograms are presented in Figures la and lb for the distance (in mm) moved by all cells in the x direction during one time unit (80 s). Negative numbers indicate a movement to the right; positive numbers indicate a movement to the left. We shall address two questions: (1) What is the probability that the granulocyte finds 0 before 2.5, given its starting position? (2) Given its starting position, what is the mean length of time until the cell reaches one of the two sides? As we shall learn, both of these questions can be easily resolved using the stochastic calculus. In addition, as Boyarsky et al. observe, the model, by its nature, suggests statistical tests. Standard Brownian motion, or the Wiener process, is a diffusion process {B’(t); t 3 O} having the following two properties: (1) E[W(t + At)- PV(t)lB’(tl= xl = 0 (2) E[(W(t + At>- W(t>>21W(t> = x1 = At
JOSEPH number
C. WATKINS
AND BIRGIT
WOESSNER
of cells
120 100 80 60 40 20
-7
-5
-3
-1
1
3
5
79tMl
1
m
(a) r lumber
of cells
160
140 120 100 80 60 40 20
-7
-5
-3
-1
3
5
79mm
(b) FIG. 1. Histograms (b) for chemotaxis.
of number
of cells versus distance
moved (a) for the control
case;
DIFFUSION
MODELS
FOR CHEMOTAXIS
279
TABLE 1 Probabilities in Control Cases of Reaching 0 before 25 Starting at x E [0,2.5] (in cm) After the Experiment Had Been Run for 96 min, and Theoretical Values for f’(a, x) for a = - 0.11 and a = 0 x 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Observed 0.964 0.930 0.896 0.863 0.828 0.791 0.751 0.707 0.663 0.621 0.579 0.535 0.490 0.445 0.403 0.363 0.324 0.287 0.249 0.212 0.173 0.132 0.090 0.045 0.000
P(-0.11,X) 0.992 0.983 0.973 0.962 0.950 0.921 0.921 0.904 0.884 0.863 0.839 0.813 0.783 0.750 0.713 0.671 0.625 0.574 0.516 0.452 0.380 0.300 0.211 0.111 0.000
P(O, x) 0.96 0.92 0.88 0.84 0.80 0.76 0.72 0.68 0.64 0.60 0.56 0.52 0.48 0.44 0.40 0.36 0.32 0.28 0.24 0.20 0.16 0.12 0.08 0.04 0.00
In words, the mean position of the Wiener process does not change in the interval of time between t and t + At, and the variance in the position is equal to At, the length of this interval. Neither of these statements is affected by the position of the process at time t. We do not expect that a cell will assume the same passive role in its environment as was taken by a pollen grain (See Brown [6] or Nelson [20]). The cell may have a preferred direction, particularly in the presence of an attractant, and this preference may depend upon its position. The variance of the motion may also depend on the cell’s position if the nature of its activity depends on its position. Letting p(x) and U’(X) be functions that represent this preference in direction and level of activity, we can again use conditional expectation to quantify cell movement. Let X(t) be the position of the cell at time t, and let At be a short interval of time. The discussion
280
JOSEPH
C. WATKINS
TABLE
AND BIRGIT
WOESSNER
2
Probabilities in Chemotaxis of Reaching 0 before 25 Starting at x E [0,25] (in cm) After the Experiment Had Been Run 80 min, and Theoretical Values for p(a, x) for a = -3.15 and a = -0.12 X
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Observed
PC-3.15,X)
P(-0.12,X)
0.999
0.993 0.986 0.977 0.968 0.957 0.945 0.931 0.912 0.898 0.878 0.856 0.831 0.803 0.771 0.735 0.695 0.649 0.598 0.540 0.475 0.401
0.992 0.984 0.975 0.965 0.955 0.944 0.932 0.917 0.901 0.884 0.883 0.840 0.813 0.783 0.746 0.703 0.655 0.601 0.541 0.472 0.391 0.304 0.209 0.110 0.000
above can be summarized
0.998 0.957 0.000
0.318 0.225 0.119 0.000
by stating that
and that
In other words, over a short period of time, the mean and variance of the displacement are nearly multiples of the length of the time interval. As the cell moves away from x, the function values for p and (T may change, and the approximations become worse. The o(At) indicates that a correction is necessary but that the size of the correction has an order smaller than At,
DIFFUSION
MODELS
281
FOR CHEMOTAXIS
the length of this time interval. Therefore, the motion at the point x can be closely approximated by Brownian motion. The parameters are not 0 and 1, as in the case of the Wiener process, but rather p(x) and U*(X). Thus, to estimate p and (T, fix 5 E[O, L] as the starting point for cell movement. Take n measurements of the horizontal cell displacements of n cells starting at 5. Let xi, x2,.. ., X, be the x coordinate at a time At later. Then the standard estimator for drift coefficient /L is
and for the diffusion
coefficient
G2(5)= j& In the Boyarski-Noble-Peterson Note that
E _(&) [
f12 the standard
estimator
is
k [~~-5-fi(~)At]~. i=l
experiment,
L = 25 cm and At = 80 s.
[X(t+At)-X(t)-p(X(t))At]lX(t)=*l_o(At)
and
E
gt;ct))[X(t+At)-X(t)-~(X(t))At]
)‘lX(t)=x 1
= At + o(At), and the statistics for the short time displacement
+&)) [X(t+At)-X(t)-~L(X(t))Atl are nearly the same as those for the Wiener $(t))
process. Therefore,
[X(t+At)-X(t)-p(X(t))At]=W(t+At)-W(t)
and X(t+At)-X(t)=p(X(t))At+v(X(t))[W(t+At)-W(t)].
282
JOSEPH C. WATKINS
In the modeling of ideas justify the use Differential calculus change in the motion equation above using
AND BIRGIT
WOESSNER
systems having no randomness, similar approximation of a differential equation to model the motion. is involved because we are dealing with rates of as At -+ 0. Here we shall do the same by writing this stochastic differentials:
dX(t)=p(X(t))dt
+c(X(t))dW(t).
We shall call this equation a stochastic differential equation. The function CLis called the drift coefficient, and u2 is called the diffusion coefficient. In their initial suggestion that diffusions be used to model chemotaxis, Keller and Segel[17] use the terms motility for (T and chemotaxis coefficient for p. Properties of the solution can be determined by finding the mean value of properly chosen functions of the solution. In order to perform these operations, we must develop a calculus for these stochastic integrals. First let’s review the change-of-variables formula for ordinary integrals. If f and g are smooth functions, then the usual change-of-variables formula for integrals uses the chain rule,
f( g( t)) = f(s(O>) + or, in differential
$f’(u(4)dds)
form,
df(g(t)) =f’(g(t))ddt). If g is replaced by the Wiener process, then this formula no longer holds. The fluctuation of the Wiener process is too rapid to stop with one term in this Taylor’s series expansion. We must include the second term.
=f’(W(t))dW(t)+~f”(W(t))[dW(t)l’.
df(Wt)) Because the variance
of IV(t) grows linearly,
[dW(t)l* = dt, and we have
df(W(t))=f’(W(t))dW(t)+$‘(W(t))dt.
The same formula
holds for X(t):
df( X(t)) = f’(X(t>) H(t)
+ if”“(X(t>>[~(t)l*.
To simplify this expression, we must remember to keep all differentials of order 1 or less in dt and discard all others. Because [dW(t)]* = dt, we say that the order of the differential for the Wiener process in terms of dt is
DIFFUSION MODELS FOR CHEMOTAXIS
l/2.
The following chart will remind 1
283
us of this multiplication &V(t)
dt
dt 0
dW(t) dt
rule:
0 0
Thus, [dX(t)]2={b(X(t))dt+cT(X(t))dW(t)}2 =b2(X(t))dt+2b(X(t))cr(X(t))dtdW(t) +
~2(wwv)12
=O+O+a2(W(t))dt=a2(X(t))dt.
Therefore,
we have the It6 formula
@(X(t))
= [G(t))f’(X(t)) + @(x(t)) = (U+
+ +~2(X(t))frr(X(t))]
dt
dW(t)
ta2Y’)( x(t))
dt + u(f’)(
x(t))
dW(t)),
where the notation of the second line means to multiply the functions together and evaluate them at X(t). We now write this expression in integral form to obtain
f(x(t))
+(x(o))
=j,l(pf’+~~*f”)(X(s)) d~+lnl(uf’)(X(‘))dW(9).
We shall call integration with respect to the Wiener process an It8 integral. The stochastic calculus can accommodate a certain class of random times called stopping times or Markov times. A random time T is called a Markov time for the stochastic process X(t) if its occurrence can be determined without looking at the process beyond the time of occurrence. Therefore, the time when a granulocyte exits a region can be determined by watching its movement up to this time. We need not continue watching. However, the last time a granulocyte cell exits a region cannot be determined with certainty unless we continue watching. Such random times are not Markov times, because they anticipate the future. Now, we apply the It8 formula to the stochastic integral equation with a Markov time.
284
JOSEPH C. WATKINS
AND BIRGIT
WOESSNER
It8 integrals have expected value zero. This fact holds also in the case where the interval of integration ends in a Markov time. Thus,
E[f(XW) -f(X(O))IX(O) = xl
&f’+
;g2f”)( X(s))
Generally, we shall write the conditional abbreviated form E’, and thus we have
dslX(0)
expectation
=
x].
E[. IX(O) = x] in the
This is our basic equation and is known as Dynkin’s formula. It applies whenever the function f is bounded and has two continuous derivatives and the Markov time T has a finite mean. (See Arnold [2]). Turning to question 1, let P(x) denote the probability that a cell starting at x E [0, L] reaches zero before it reaches L, and let f be the solution to the differential equation
with boundary
conditions and
f(0) = 1
f(L)=O.
Let T denote the random time when the cell first reaches zero or L. Then r is a Markov time, and Dynkin’s formula applies.
=
E”
f(x)+
[
j,;Ods
1
=f(x) because f was chosen to solve the differential equation. Note that X(T) can take only two values, 0 and L, depending upon which end the cell first reaches. Thus,
f (X(T))
=
1
ifX(7)=0
i 0
ifX(7)=L
DIFFUSION
285
MODELS FOR CHEMOTAXIS
by the boundary
conditions.
Therefore,
E”[f(X(7))]=1F{f(X(7))=1}+0F(f(X(7))=0} = P”{ X( 7) = 0} = P”( cell reaches 0 before L) = P(x). Here Px denotes Thus,
the probability
f(x)
=
conditioned
-wfw~))1=
on the event that X(O)= x.
P(x).
In other words, P is the solution to the differential equation above. If uz is positive and continuous, then the solution to the differential equation, which is easy to verify directly, is
with =exp
S(y)
[
-/r:sd[].
Now an estimate based on the model can be made for P. We just replace the integral used to define S with a Riemann sum for this integral using the estimators for the drift and the diffusion. Thus, if we estimate k and (T* at the points t,, &, . . . , tN, with 0 = to < 5, < . . . < lN = L, then we take
for y E (sj, s,,,].
The estimator
for P will be
N-l
N-l
C
g(Sj)(Sj+l-
j-k
tj>
C
s^(Cj)(tj+*-Sj)
j=O
In the model under study, ti = i cm and N = 25. In the neutral cases, we anticipate that p(x) = 0 for all x. In this instance, P is the linear function P(x)=(L-x)/L.
286
JOSEPH
C. WATKINS
AND BIRGIT
WOESSNER
1.0
0.6
x
Observed
0.6 b g 1 e P
0.4
0.2
0.0 20
10
0
Centimeters
FIG. 2. Observed
vs. theoretical
probabilities
under
control
conditions.
Boyarsky et al. [5] state that the nearness of the observed probabilities to linearity is evidence against chemotaxis in the control case and indeed, as Figure 2 shows us, their observations are very close to linear. Our preferred line of reasoning is to take local observations and a model to make global predictions and base our reliability on the model on the accuracy of these predictions. See Woessner [35]. Nevertheless, their line of argument is useful, and we shall explore it. Assume for the moment that 2~/(+ is some constant but unknown parameter a. In this case, the hitting probabilities take the simpler form
We shall use a weighted
F(a)=
where
n,
is the number
least squares estimate
for a by minimizing
[fYa>x)-P(n)]z *=,n,p(u,x)[l-P(a,x)l’ :
of cells originating
in strip x and P(X) is the
DIFFUSION MODELS FOR CHEMOTAXIS
287
observed probability of hitting weights P(a, x)[l- P(a, x)1/n,
the
zero starting at x. We have chosen for the following reason. Define
if the ith cell in strip x reaches 0 before L if the ith cell in strip x reaches L before 0. The observations {fii(x); 1~ i < n,, 1~ x < 24) are assumed to be independent Bernoulli random variables with parameter P(a, x). Then
Under this assumption p(x) has mean P(a,x) and variance P(u,x)[lP(u,x)]/n,. Thus the weights standardize the variance at 1. With 90-120 cells in each experiment and with the experiment repeated six times, we have n x = 30. Using the data in Table 1, we find that a = -0.0066 minimizes F(u) to have the value 0.36. We can also use the weighted least squares fit as a x2 test for the null hypothesis He:
the parameter
a = 0.
The x2 statistic is much too small to reject this hypothesis. If we apply these ideas to the experiments where bacteria are present, then we find that a = 0.12 minimizes F(u) at the value 0.25. F(u) has a x2 distribution with 24 degrees of freedom, and x&0.95) = 36.4. As Figure 2 shows, the fit of the data to this model is strikingly good. A 95% confidence interval for a is [ -0.185, - 0.061, meaning that each point in this interval has a x2 value less than 36.4. Also, the x2 test rejects Ha. Something beyond pure constant random mobility is occurring. Unfortunately, problems with the data prohibit us from using these experiments to give a complete application of our methodology. To describe quickly the source of the difficulties, focus on the histogram in Figure la for the control case. Let d,, d,,.. . be the change in the x coordinate for each of the individual cells. Assume that these displacements are independent and identically distributed. Let p be the mean of their displacements and C* the variance. Then
2 d,-np z,=
i=l
(nay2 has an approximately
normal
distribution
with mean 0, variance
1. In this
288
JOSEPH C. WATKINS AND BIRGIT WOESSNER
0.6 -
X
Observed
0.0 I 0
10
20
Centimeters
FIG. 3. Observed vs. theoretical probabilities for chemotaxis.
case, we take as our null hypothesis H,:
p = 0.
All displacements smaller than 3 mm have been considered zero, and thus we consider, for the purpose of testing the hypothesis, only those displacements that are 3 mm or greater. We obtain n = 144
and
2 di= -208. i=l
We take the sample variance I?* = 21.55 for the unknown variance u2. Substituting these values, we find that Z, = -3.73. A normal random variable with mean 0 and variance 1 takes on a value less than or equal to -3.09 with probability 0.001. On this basis, we reject H,. The data presented are unlikely to be a sample from a distribution having mean 0. These data give 2&/G2 = -0.11 in the control case, and for chemotaxis -3.15. The predicted probabilities for these parameter values 2fi/62= are tabulated in Tables 1 and 2 and graphed in Figures 3 and 4. A second source of concern is seen in the histogram for chemotaxis (Figure lb). In
DIFFUSION
MODELS
FOR
289
CHEMOTAXIS
1.c
P(-3.15.X) 06
0.6
F
:Fj
2
e
0.
0.4
0.2
0.0 10
0
20
Centimeters FIG. 4. Estimates
for n obtained
from histogram.
this distribution, we find large displacements away from the attractant, displacements that we did not find in the control environment. Despite these problems, we see that taking 2p/u2 constant gives excellent agreement with the data. Gruler and Biiltman 1141give an experimental value for & 2: 240 pm*/min and see a persistence time on the order of a minute. In the Boyarsky-Noble-Peterson experiment we see nearly normal distributed data after 80 s and consequently no evidence that persistence lasts longer than 80 s. Persistence times this short are unusual and may be an artifact of restricting displacements under consideration to those that are at least 3 mm. (We thank a referee for this last remark.) As we shall see in the discussion of approximation theorems, 2p/a2 can be constant without p or CTbeing constant. In order to separate this ratio, we now turn to calculating the mean length of time required for the granulocyte to reach one of the two sides. Let T(x) denote the expected time that a cell starting at x E 10, Ll reaches 0 or L, and let f be the solution to the differential equation /.L(x)f’(
x) +
+a*( x)f”(X)
= - 1
290 with boundary
JOSEPH C. WATKINS
AND BIRGIT
WOESSNER
conditions f(0) = 0
f(L) =o.
and
Let 7, as before, be the random time when the cell first reaches Again, we write Dynkin’s formula, ’
0 or L.
= f(x)+EQp)dr] =f(x)-E”[ps] = f(x)-
EX7.
As before, X(T) can take only two values, 0 and L. In both instances, f(X(T)) = 0 by the boundary conditions. Therefore, 0= f(X)-EXT or
Thus, T(x) solves the differential then
equation
stated above. Let g(x) = f’(x>;
j_+)g(x)+&+)g’(x)
= -1
is a first-order linear nonhomogeneous ordinary differential a’(x) is continuous and positive, then we divide u* to obtain
The integrating
factor for this differential S(y)=exp
equation
If
is
-~y$!$fdt]. [
Thus, we integrate. The constants of integration boundary conditions, and we obtain
W=@s)
equation.
can be found
dZ)dY, ic-joYu2(z;s(z)
using the
DIFFUSION
MODELS
291
FOR CHEMOTAXIS
where the constant -1
c= i ~L~(Y)~y~2(z;s(z)
dZdY)(kLSCY)dY)
.
Thus, the model and the local estimates for I_Cand a2 at the points O=&)
Unicellular Movement
JOSEPH C. WATKINS Department of Mathematics, University of Southern California, Los Angeles, California 90089-1113 AND BIRGIT WOESSNER Fachbereich Mathematik, Universitiit Urn, Federal Republic of Germany Received 7 September 1989; revised 30 November 1990
ABSTRACT A program is developed for applying stochastic differential equations to models for chemotaxis. First a few of the experimental and theoretical models for chemotaxis both for swimming bacteria and for cells migrating along a substrate are reviewed. In physical and biological models of deterministic systems, finite difference equations are often replaced by a limiting differential equation in order to take advantage of the ease in the use of calculus. A similar but more intricate methodology is developed here for stochastic models for chemotaxis. This exposition is possible because recent work in probability theory gives ease in the use of the stochastic calculus for diffusions and broad applicability in the convergence of stochastic difference equations to a stochastic differential equation. Stochastic differential equations suggest useful data for the model and provide statistical tests. We begin with phenomenological considerations as we analyze a onedimensional model proposed by Boyarsky, Noble, and Peterson in their study of human granulocytes. In this context, a theoretical model consists in identifying which diffusion best approximates a model for cell movement based upon theoretical considerations of cell phsyiology. Such a diffusion approximation theorem is presented along with discussion of the relationship between autocovariance and persistence. Both the stochastic calculus and the diffusion approximation theorem are described in one dimension. Finally, these tools are extended to multidimensional models and applied to a threedimensional experimental setup of spherical symmetry.
1.
INTRODUCTION
In 1905, Einstein [ll] gave his explanation for Brownian motion via a kinetic theory of molecules in part to convince the vitalists once and for all of the hopelessness of their point of view. In a strange twist of irony scarcely a decade after Einstein’s description of the subject, Przibaum [27] MATHEMATICAL
BIOSCIENCES
OElsevier Science Publishing 655 Avenue of the Americas,
104:271-303
(1991)
Co., Inc., 1991 New York, NY 10010
271 0025-5564/91/$03.50
272
JOSEPH C. WATKINS
AND BIRGIT
WOESSNER
found that Brownian movement was an adequate description for the movement of unicellular organisms. (See Thompson [31], page 76). Seventy years later, biologists have much more technologically advanced tools than a microscope, a pencil, and a metronome to investigate the motion of such organisms. At the same time, mathematicians have much more technologically advanced tools to understand Brownian motion and related random processes. The purpose of this paper is to inject some of these probabilistic tools into the biological literature in such a way that biologists can avail themselves of them. In conformity with recent mathematical models, we consider two types of cellular movement: the swimming motion of bacteria and the migrating motion of cells moving along a substrate. These models can be used to try to understand cell movement from the more superficial phenomenological point of view. In other words, the experimenter can take data and construct a consistent model with no real regard for the underlying cellular or molecular processes that might give rise to the activities under study. More fundamentally, we can propose that the randomness in the mathematical model derives from intrinsic fluctuations in the local environment and in the sensing mechanism. Cell movement is thus modeled as a stochastic process whose randomness derives from these two sources. In the particular case of Escherichiu coli, a detailed statistical description of the movement was presented by Berg and Brown [3]. By building a microscope that automatically follows the individual cells, they found an alternating sequence of two types of behavior quite different in character. In the first the change in direction is gradual, and in the second the change in direction is abrupt. They called these two behavior types runs and twiddles. During a run the flagellar filaments of a bacterium work together in a bundle as it moves. When the cell begins to twiddle, the bundle loosens. The bundle then re-forms, and the cell travels in a new direction that depends upon the change in orientation of the bundle relative to the body of the cell. Berg and Brown observed that the change in direction between the run preceding and the one following a twiddle is an angle whose distribution does not vary from twiddle to twiddle. Also, the speed of the bacterium is essentially constant during a run and between two different runs. Twiddle times are short compared to run times, and in an isotropic environment both the twiddle times and the run times have an exponential distribution. If an attractant is introduced into the environment, the distribution of twiddle times remains unchanged. The run times still have an exponential distribution. However, the parameter in the distribution depends on the angle between the run and the gradient of the concentration of the attractant. The bacterium seems to have no control during a twiddle. Consequently, if a cell has a preferred direction, it can influence its movement only by altering the length of the run times.
DIFFUSION
MODELS
FOR CHEMOTAXIS
273
Polymorphonuclear leukocyte or granulocyte cells will serve as an example of a cell whose typical behavior is locomotion along a surface. In actuality, these cells thrust themselves forward by means of a lamellipodium. The forming of each lamellipodium involves the streaming of cytoplasm, which brings the bulk of the cell forward. Trinkhaus [33] describes this mode of locomotion in great detail and speculates on the mechanisms that might explain its behavior. A comprehensive statistical description of leukocyte movement will be quite difficult. Under optimal conditions, Zigmond et al. [36] describe the leukocyte as morphologically and behaviorally polarized. Movement is smooth and is governed by a stable lamellipodium. Under suboptimal or artificial conditions, leukocytes can make erratic changes in polarity due to competition among multiple lamellipodia (see Gerisch and Keller 1131 and Shields and Haston [291X How does one turn these observations into a mathematical model for chemotaxis? As usual, one wants to simplify the situation to produce a model that captures the essential features of the actual phenomena in a way that can be treated mathematically. Ideally, the model should be susceptible to statistical tests that check its reliability. To this date, most models of swimming bacteria assume that a run is linear motion at a constant speed and a twiddle is a change in direction. In addition, models typically assume that all cells have the same pattern of behavior and that the cells do not interact, for example via contact or chemical transmitters. When experimenters try to conduct experiments in an environment that is as homogeneous (in time or in space) as possible, the model reflects this attempt by assuming a suitably homogeneous environment. More detailed models include the exponential distribution for the random lengths of twiddle times and run times. The usual models following from these assumptions are random walks, or more generally Markov processes. The Berg-Brown experiment is prototypical for unicellular movement. Data are collected for short-time behavior and are synthesized to supply the parameters in the model. The statistical checks for the reliability of the model involve its ability to predict long-term behavior. Typical questions might include: Do the cells have a preferred direction in the presence or the absence of an attractant? What is the probability that the cell finds the attractant? How much time, on the average, does the cell require to find the attractant? Peterson and Noble [26] use a two-dimensional random walk model to analyze human granulocyte movement. If this random walk has no preferred direction-interpreted to mean that the direction of the run has mean zero-then the variance in the position of a cell should be linearly proportional to the number of runs. Peterson and Noble observed the movement of granulocytes in the absence of bacteria and subsequently in the presence of bacteria. Their data were the number of runs and the
274
JOSEPH C. WATKINS
AND BIRGIT
WOESSNER
displacement during each of the runs. The total distance traveled was computed, and the results were analyzed by performing a t-test to determine whether the mean direction of a run was indeed zero. In the absence of bacteria, the difference between the actual distance and the predicted distance was insignificant (P > O.OS>,whereas when bacteria were present, the difference was highly significant (P < 0.001). Their conclusion: The runs have mean displacement zero if no bacteria are present. However, if bacteria are present the data are extremely unlikely to be a sample of random displacements having mean zero. They interpret this unlikeliness as evidence for chemotaxis. Nossal and Weiss [21] propose a descriptive theory for cell movement in two dimensions. Assuming that the starting point for the cell is the origin of coordinates and that the attractant is located on the x axis, they derive expressions for the mean and variance of the cell displacement after n runs. Their model has considerable detail: Run and twiddle times are exponential in the manner described by Berg and Brown, the turning angle is uniformly distributed, and path lengths and turning angles are independent. Using the words of Berg and Brown, the model has an abundance of Poisson statistics, and consequently the model is essentially Markovian. A random process is called Markovian if the probability distribution of any future occurrence depends on the history of the process up to the present time only via the present situation. The exponential distribution is a typical feature in Markov processes because the probability per unit time of the termination of a twiddle or the termination of a run is constant. The past, in this case the length of the present run or twiddle, is no help in predicting when the run or twiddle will end. Thus the pair (X(t),V(t)) taken together form a Markov process. If we see the cell at X(t) moving with velocity V(t), we know that it will continue in this direction for an exponential length of time and that the parameter for this exponential random variable depends on X(t) and V(t). As we begin to make predictions about the future after time t, we see that we need only X(t) and V(t) to make these predictions. To date, Wolfgang Alt [l] has developed the most comprehensive axiomatic model for chemotaxis. His model begins with a general setup for the two-parameter stochastic process (X(t), V(t)) as a time inhomogeneous semi-Markov process. As we have seen in a Markov process, V(t) is constant for an exponential length of time. The parameter in this exponential distribution can depend on position and on time. A semi-Markov process maintains the Markov property for the sequence of velocities v,, vz, v,, . . . that a cell happens to choose, but allows the run time to take on an arbitrary distribution. As before, these run times are independent and have an identical distribution, which can depend on the location of the cell and on time. (See Othmer et al. [24] for more on this. They call this
DIFFUSION MODELS FOR CHEMOTAXIS
275
process a kangaroo process or a position jump process. See Cinlar [81 for an extensive development of semi-Markov processes.) Alt soon restricts his attention to time-homogeneousprocesses. This acknowledges the attempts that experimenters make to create a time-homogeneous environment and the difficulties in statistical inference for time-inhomogeneous processes but ignores the fact that cells fatigue or age. From biochemical considerations, Ah can further restrict his model to give plausible expressions for the sequence of velocities and for the times for which the cell maintains these velocities. For more, see Okubo [23]. In this paper, we adopt a point of view that is part of the modern trend in probabilistic analysis; that is, we focus on the random processes themselves rather than on the partial differential operators, difference equations, or integral kernels that are featured in more analytical treatments. The hope is that by taking this point of view we can find a more unified framework for the present discussion and a more natural path for generalization in future discussions. After all, the objects of study are random processes and thus they rightfully belong at the center of attention. Let us emphasize that this treatment supplements, but cannot replace, the sound analytical treatment that, for example, Wolfgang Alt has given. Indeed, Alt’s analysis is essential in identifying parameters in the model and in estimating the errors involved in using a diffusion approximation. The model advocated in this manuscript is a diffusion model. The questions that a biologist should ponder before adopting a diffusion model are two:’ (1) Is the motion continuous? (2) Does the distribution of future motion of the cell depend only on its present position? Typically, we can easily decide on an answer to the first question. An affirmative answer to question 2 is always rendered with some reservations. For example, we have just argued that direction is also needed to predict future movement in the Nossal-Weiss model. However, one may reason that bacteria make many runs, with a direction change at the end of each run. So the current direction, which the cell follows for only a very short time, has very little impact on the cell’s ability to find its attractant. If this rationale is convincing, then, for the purpose of modeling, yes is a reasonable response to question 2. Certainly, in the model described above, this response is mathematically well-grounded. We are asking a question concerning the result of a large number of small random fluctuations, and in such an instance the limit theorems of probability theory play a key role. Answering yes for question 2 is similar to accepting the sample mean of random variables as the mean and then ascribing a normal distribution about the mean. In other words, a
276 yes
JOSEPH C. WATKINS AND BIRGIT WOESSNER
to question 2 reflects the fact that the law of large numbers and the central limit theorem can be used to give good approximations. In the mathematical literature, the Nossal-Weiss model is called to transport process (see Chandresekar [7]), and answering the two questions in the affirmative would constitute an affirmation of the applicability of a version of the central limit theorem appropriate in this context. Because the details of the velocity term vanish in the limit and are replaced by average values for the velocity, the limit theorem is sometimes called homogenization. In a homogeneous environment this average quantity appears as the slope of a linearly increasing variance. In following the tracks of paramecia, Przibaum witnessed this behavior “with great exactitude.” On the other hand, this linearity was not found with great exactitude by Gruler and Biiltman 1141in their study of leukocytes. However, linear growth in the variance of displacement occurred after a period of sublinear growth. This is one of many aspects of a behavior properly called persistence. As we shall see, if the ratio between the typical time of events of interest and this persistence time is sufficiently large, then a diffusion model is appropriate. Consequently, if the variance is not asymptotically linear, then a diffusion model is never appropriate. For example, if the variance grows at some power of time, we may choose to investigate a fractal model. The definition of a diffusion is a very simple statement: A diffusion is a Markov process whose motion is continuous. Not all diffusions can be easily estimated statistically. However, if we can closely approximate the motion from each point by a Brownian motion that may not drift in some direction, then, as we shall learn, the diffusion process may be represented as the solution to a stochastic differential equation. Chemotaxis, the physical process of interest, is a multidimensional process. The theory we shall describe extends easily to higher dimensions and is complicated only by the additional notation necessary in multidimensional calculus. Although the diffusion equation appears frequently in the literature on chemotaxis, there has been, to date, no systematic exploitation of the probabilistic perspective. Here, we formulate the problem in terms of a stochastic differential equation. This formulation is neither more nor less than the classical diffusion equation but serves to introduce us to a new set of mathematical tools. For example, Boyarsky et al. [5] mention the use of stochastic differential equations. However, in the absence of a probabilistic description of the events, they cannot devise good statistical tests to evaluate their data on boundary-crossing probabilities. These boundarycrossing probabilities provide only one of several possible statistical and experimental methods upon which one may decide on a diffusion model. For example, the mean time until this boundary crossing is easily computed and provides more evidence for accepting or rejecting any given model. If the experiment can be run in several geometries, then we can be more
DIFFUSION
MODELS
FOR CHEMOTAXIS
277
certain that our data are revealing something intrinsic concerning cell movement rather than some artifact of the geometry. Thus in Section 4 we give an example with spherical symmetry to demonstrate the adaptability of the method. With this case study, it is hoped that the experimenter can establish the appropriate stochastic differential equation and the corresponding statistical tests. 2.
ANALYSIS
OF A ONE-DIMENSIONAL
MODEL
We apply these ideas to the situation in which a stochastic differential equation model was first proposed. In [S], Boyarsky et al. observed the motion of human granulocytes. Ninety to 120 cells were placed in a capillary tube, and their movement was recorded by photographing a microscopic image (54 X) every 8 s. The x component of the displacement was recorded at 80 s, the shortest time in which movement could be reliably measured. Each photograph was divided into 2.5 strips 1 cm wide, and all cells were identified with their original strip. This experiment was reproduced six times in a neutral case (control) and six times with bacteria present on the left side of the capillary tube (chemotaxis). These data are summarized in Figure 1 with histograms. A move of less than one-half the cell diameter (3 mm under magnification) is difficult to detect and hence is set to zero on the histogram. A second set of data recorded the probability of reaching zero, the left side, before 2.5, the right side of the photograph at later times. The data recorded after an experiment had run 96 min in chemotaxis and 80 min in the control case are included in Tables 1 and 2, respectively. Histograms are presented in Figures la and lb for the distance (in mm) moved by all cells in the x direction during one time unit (80 s). Negative numbers indicate a movement to the right; positive numbers indicate a movement to the left. We shall address two questions: (1) What is the probability that the granulocyte finds 0 before 2.5, given its starting position? (2) Given its starting position, what is the mean length of time until the cell reaches one of the two sides? As we shall learn, both of these questions can be easily resolved using the stochastic calculus. In addition, as Boyarsky et al. observe, the model, by its nature, suggests statistical tests. Standard Brownian motion, or the Wiener process, is a diffusion process {B’(t); t 3 O} having the following two properties: (1) E[W(t + At)- PV(t)lB’(tl= xl = 0 (2) E[(W(t + At>- W(t>>21W(t> = x1 = At
JOSEPH number
C. WATKINS
AND BIRGIT
WOESSNER
of cells
120 100 80 60 40 20
-7
-5
-3
-1
1
3
5
79tMl
1
m
(a) r lumber
of cells
160
140 120 100 80 60 40 20
-7
-5
-3
-1
3
5
79mm
(b) FIG. 1. Histograms (b) for chemotaxis.
of number
of cells versus distance
moved (a) for the control
case;
DIFFUSION
MODELS
FOR CHEMOTAXIS
279
TABLE 1 Probabilities in Control Cases of Reaching 0 before 25 Starting at x E [0,2.5] (in cm) After the Experiment Had Been Run for 96 min, and Theoretical Values for f’(a, x) for a = - 0.11 and a = 0 x 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Observed 0.964 0.930 0.896 0.863 0.828 0.791 0.751 0.707 0.663 0.621 0.579 0.535 0.490 0.445 0.403 0.363 0.324 0.287 0.249 0.212 0.173 0.132 0.090 0.045 0.000
P(-0.11,X) 0.992 0.983 0.973 0.962 0.950 0.921 0.921 0.904 0.884 0.863 0.839 0.813 0.783 0.750 0.713 0.671 0.625 0.574 0.516 0.452 0.380 0.300 0.211 0.111 0.000
P(O, x) 0.96 0.92 0.88 0.84 0.80 0.76 0.72 0.68 0.64 0.60 0.56 0.52 0.48 0.44 0.40 0.36 0.32 0.28 0.24 0.20 0.16 0.12 0.08 0.04 0.00
In words, the mean position of the Wiener process does not change in the interval of time between t and t + At, and the variance in the position is equal to At, the length of this interval. Neither of these statements is affected by the position of the process at time t. We do not expect that a cell will assume the same passive role in its environment as was taken by a pollen grain (See Brown [6] or Nelson [20]). The cell may have a preferred direction, particularly in the presence of an attractant, and this preference may depend upon its position. The variance of the motion may also depend on the cell’s position if the nature of its activity depends on its position. Letting p(x) and U’(X) be functions that represent this preference in direction and level of activity, we can again use conditional expectation to quantify cell movement. Let X(t) be the position of the cell at time t, and let At be a short interval of time. The discussion
280
JOSEPH
C. WATKINS
TABLE
AND BIRGIT
WOESSNER
2
Probabilities in Chemotaxis of Reaching 0 before 25 Starting at x E [0,25] (in cm) After the Experiment Had Been Run 80 min, and Theoretical Values for p(a, x) for a = -3.15 and a = -0.12 X
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Observed
PC-3.15,X)
P(-0.12,X)
0.999
0.993 0.986 0.977 0.968 0.957 0.945 0.931 0.912 0.898 0.878 0.856 0.831 0.803 0.771 0.735 0.695 0.649 0.598 0.540 0.475 0.401
0.992 0.984 0.975 0.965 0.955 0.944 0.932 0.917 0.901 0.884 0.883 0.840 0.813 0.783 0.746 0.703 0.655 0.601 0.541 0.472 0.391 0.304 0.209 0.110 0.000
above can be summarized
0.998 0.957 0.000
0.318 0.225 0.119 0.000
by stating that
and that
In other words, over a short period of time, the mean and variance of the displacement are nearly multiples of the length of the time interval. As the cell moves away from x, the function values for p and (T may change, and the approximations become worse. The o(At) indicates that a correction is necessary but that the size of the correction has an order smaller than At,
DIFFUSION
MODELS
281
FOR CHEMOTAXIS
the length of this time interval. Therefore, the motion at the point x can be closely approximated by Brownian motion. The parameters are not 0 and 1, as in the case of the Wiener process, but rather p(x) and U*(X). Thus, to estimate p and (T, fix 5 E[O, L] as the starting point for cell movement. Take n measurements of the horizontal cell displacements of n cells starting at 5. Let xi, x2,.. ., X, be the x coordinate at a time At later. Then the standard estimator for drift coefficient /L is
and for the diffusion
coefficient
G2(5)= j& In the Boyarski-Noble-Peterson Note that
E _(&) [
f12 the standard
estimator
is
k [~~-5-fi(~)At]~. i=l
experiment,
L = 25 cm and At = 80 s.
[X(t+At)-X(t)-p(X(t))At]lX(t)=*l_o(At)
and
E
gt;ct))[X(t+At)-X(t)-~(X(t))At]
)‘lX(t)=x 1
= At + o(At), and the statistics for the short time displacement
+&)) [X(t+At)-X(t)-~L(X(t))Atl are nearly the same as those for the Wiener $(t))
process. Therefore,
[X(t+At)-X(t)-p(X(t))At]=W(t+At)-W(t)
and X(t+At)-X(t)=p(X(t))At+v(X(t))[W(t+At)-W(t)].
282
JOSEPH C. WATKINS
In the modeling of ideas justify the use Differential calculus change in the motion equation above using
AND BIRGIT
WOESSNER
systems having no randomness, similar approximation of a differential equation to model the motion. is involved because we are dealing with rates of as At -+ 0. Here we shall do the same by writing this stochastic differentials:
dX(t)=p(X(t))dt
+c(X(t))dW(t).
We shall call this equation a stochastic differential equation. The function CLis called the drift coefficient, and u2 is called the diffusion coefficient. In their initial suggestion that diffusions be used to model chemotaxis, Keller and Segel[17] use the terms motility for (T and chemotaxis coefficient for p. Properties of the solution can be determined by finding the mean value of properly chosen functions of the solution. In order to perform these operations, we must develop a calculus for these stochastic integrals. First let’s review the change-of-variables formula for ordinary integrals. If f and g are smooth functions, then the usual change-of-variables formula for integrals uses the chain rule,
f( g( t)) = f(s(O>) + or, in differential
$f’(u(4)dds)
form,
df(g(t)) =f’(g(t))ddt). If g is replaced by the Wiener process, then this formula no longer holds. The fluctuation of the Wiener process is too rapid to stop with one term in this Taylor’s series expansion. We must include the second term.
=f’(W(t))dW(t)+~f”(W(t))[dW(t)l’.
df(Wt)) Because the variance
of IV(t) grows linearly,
[dW(t)l* = dt, and we have
df(W(t))=f’(W(t))dW(t)+$‘(W(t))dt.
The same formula
holds for X(t):
df( X(t)) = f’(X(t>) H(t)
+ if”“(X(t>>[~(t)l*.
To simplify this expression, we must remember to keep all differentials of order 1 or less in dt and discard all others. Because [dW(t)]* = dt, we say that the order of the differential for the Wiener process in terms of dt is
DIFFUSION MODELS FOR CHEMOTAXIS
l/2.
The following chart will remind 1
283
us of this multiplication &V(t)
dt
dt 0
dW(t) dt
rule:
0 0
Thus, [dX(t)]2={b(X(t))dt+cT(X(t))dW(t)}2 =b2(X(t))dt+2b(X(t))cr(X(t))dtdW(t) +
~2(wwv)12
=O+O+a2(W(t))dt=a2(X(t))dt.
Therefore,
we have the It6 formula
@(X(t))
= [G(t))f’(X(t)) + @(x(t)) = (U+
+ +~2(X(t))frr(X(t))]
dt
dW(t)
ta2Y’)( x(t))
dt + u(f’)(
x(t))
dW(t)),
where the notation of the second line means to multiply the functions together and evaluate them at X(t). We now write this expression in integral form to obtain
f(x(t))
+(x(o))
=j,l(pf’+~~*f”)(X(s)) d~+lnl(uf’)(X(‘))dW(9).
We shall call integration with respect to the Wiener process an It8 integral. The stochastic calculus can accommodate a certain class of random times called stopping times or Markov times. A random time T is called a Markov time for the stochastic process X(t) if its occurrence can be determined without looking at the process beyond the time of occurrence. Therefore, the time when a granulocyte exits a region can be determined by watching its movement up to this time. We need not continue watching. However, the last time a granulocyte cell exits a region cannot be determined with certainty unless we continue watching. Such random times are not Markov times, because they anticipate the future. Now, we apply the It8 formula to the stochastic integral equation with a Markov time.
284
JOSEPH C. WATKINS
AND BIRGIT
WOESSNER
It8 integrals have expected value zero. This fact holds also in the case where the interval of integration ends in a Markov time. Thus,
E[f(XW) -f(X(O))IX(O) = xl
&f’+
;g2f”)( X(s))
Generally, we shall write the conditional abbreviated form E’, and thus we have
dslX(0)
expectation
=
x].
E[. IX(O) = x] in the
This is our basic equation and is known as Dynkin’s formula. It applies whenever the function f is bounded and has two continuous derivatives and the Markov time T has a finite mean. (See Arnold [2]). Turning to question 1, let P(x) denote the probability that a cell starting at x E [0, L] reaches zero before it reaches L, and let f be the solution to the differential equation
with boundary
conditions and
f(0) = 1
f(L)=O.
Let T denote the random time when the cell first reaches zero or L. Then r is a Markov time, and Dynkin’s formula applies.
=
E”
f(x)+
[
j,;Ods
1
=f(x) because f was chosen to solve the differential equation. Note that X(T) can take only two values, 0 and L, depending upon which end the cell first reaches. Thus,
f (X(T))
=
1
ifX(7)=0
i 0
ifX(7)=L
DIFFUSION
285
MODELS FOR CHEMOTAXIS
by the boundary
conditions.
Therefore,
E”[f(X(7))]=1F{f(X(7))=1}+0F(f(X(7))=0} = P”{ X( 7) = 0} = P”( cell reaches 0 before L) = P(x). Here Px denotes Thus,
the probability
f(x)
=
conditioned
-wfw~))1=
on the event that X(O)= x.
P(x).
In other words, P is the solution to the differential equation above. If uz is positive and continuous, then the solution to the differential equation, which is easy to verify directly, is
with =exp
S(y)
[
-/r:sd[].
Now an estimate based on the model can be made for P. We just replace the integral used to define S with a Riemann sum for this integral using the estimators for the drift and the diffusion. Thus, if we estimate k and (T* at the points t,, &, . . . , tN, with 0 = to < 5, < . . . < lN = L, then we take
for y E (sj, s,,,].
The estimator
for P will be
N-l
N-l
C
g(Sj)(Sj+l-
j-k
tj>
C
s^(Cj)(tj+*-Sj)
j=O
In the model under study, ti = i cm and N = 25. In the neutral cases, we anticipate that p(x) = 0 for all x. In this instance, P is the linear function P(x)=(L-x)/L.
286
JOSEPH
C. WATKINS
AND BIRGIT
WOESSNER
1.0
0.6
x
Observed
0.6 b g 1 e P
0.4
0.2
0.0 20
10
0
Centimeters
FIG. 2. Observed
vs. theoretical
probabilities
under
control
conditions.
Boyarsky et al. [5] state that the nearness of the observed probabilities to linearity is evidence against chemotaxis in the control case and indeed, as Figure 2 shows us, their observations are very close to linear. Our preferred line of reasoning is to take local observations and a model to make global predictions and base our reliability on the model on the accuracy of these predictions. See Woessner [35]. Nevertheless, their line of argument is useful, and we shall explore it. Assume for the moment that 2~/(+ is some constant but unknown parameter a. In this case, the hitting probabilities take the simpler form
We shall use a weighted
F(a)=
where
n,
is the number
least squares estimate
for a by minimizing
[fYa>x)-P(n)]z *=,n,p(u,x)[l-P(a,x)l’ :
of cells originating
in strip x and P(X) is the
DIFFUSION MODELS FOR CHEMOTAXIS
287
observed probability of hitting weights P(a, x)[l- P(a, x)1/n,
the
zero starting at x. We have chosen for the following reason. Define
if the ith cell in strip x reaches 0 before L if the ith cell in strip x reaches L before 0. The observations {fii(x); 1~ i < n,, 1~ x < 24) are assumed to be independent Bernoulli random variables with parameter P(a, x). Then
Under this assumption p(x) has mean P(a,x) and variance P(u,x)[lP(u,x)]/n,. Thus the weights standardize the variance at 1. With 90-120 cells in each experiment and with the experiment repeated six times, we have n x = 30. Using the data in Table 1, we find that a = -0.0066 minimizes F(u) to have the value 0.36. We can also use the weighted least squares fit as a x2 test for the null hypothesis He:
the parameter
a = 0.
The x2 statistic is much too small to reject this hypothesis. If we apply these ideas to the experiments where bacteria are present, then we find that a = 0.12 minimizes F(u) at the value 0.25. F(u) has a x2 distribution with 24 degrees of freedom, and x&0.95) = 36.4. As Figure 2 shows, the fit of the data to this model is strikingly good. A 95% confidence interval for a is [ -0.185, - 0.061, meaning that each point in this interval has a x2 value less than 36.4. Also, the x2 test rejects Ha. Something beyond pure constant random mobility is occurring. Unfortunately, problems with the data prohibit us from using these experiments to give a complete application of our methodology. To describe quickly the source of the difficulties, focus on the histogram in Figure la for the control case. Let d,, d,,.. . be the change in the x coordinate for each of the individual cells. Assume that these displacements are independent and identically distributed. Let p be the mean of their displacements and C* the variance. Then
2 d,-np z,=
i=l
(nay2 has an approximately
normal
distribution
with mean 0, variance
1. In this
288
JOSEPH C. WATKINS AND BIRGIT WOESSNER
0.6 -
X
Observed
0.0 I 0
10
20
Centimeters
FIG. 3. Observed vs. theoretical probabilities for chemotaxis.
case, we take as our null hypothesis H,:
p = 0.
All displacements smaller than 3 mm have been considered zero, and thus we consider, for the purpose of testing the hypothesis, only those displacements that are 3 mm or greater. We obtain n = 144
and
2 di= -208. i=l
We take the sample variance I?* = 21.55 for the unknown variance u2. Substituting these values, we find that Z, = -3.73. A normal random variable with mean 0 and variance 1 takes on a value less than or equal to -3.09 with probability 0.001. On this basis, we reject H,. The data presented are unlikely to be a sample from a distribution having mean 0. These data give 2&/G2 = -0.11 in the control case, and for chemotaxis -3.15. The predicted probabilities for these parameter values 2fi/62= are tabulated in Tables 1 and 2 and graphed in Figures 3 and 4. A second source of concern is seen in the histogram for chemotaxis (Figure lb). In
DIFFUSION
MODELS
FOR
289
CHEMOTAXIS
1.c
P(-3.15.X) 06
0.6
F
:Fj
2
e
0.
0.4
0.2
0.0 10
0
20
Centimeters FIG. 4. Estimates
for n obtained
from histogram.
this distribution, we find large displacements away from the attractant, displacements that we did not find in the control environment. Despite these problems, we see that taking 2p/u2 constant gives excellent agreement with the data. Gruler and Biiltman 1141give an experimental value for & 2: 240 pm*/min and see a persistence time on the order of a minute. In the Boyarsky-Noble-Peterson experiment we see nearly normal distributed data after 80 s and consequently no evidence that persistence lasts longer than 80 s. Persistence times this short are unusual and may be an artifact of restricting displacements under consideration to those that are at least 3 mm. (We thank a referee for this last remark.) As we shall see in the discussion of approximation theorems, 2p/a2 can be constant without p or CTbeing constant. In order to separate this ratio, we now turn to calculating the mean length of time required for the granulocyte to reach one of the two sides. Let T(x) denote the expected time that a cell starting at x E 10, Ll reaches 0 or L, and let f be the solution to the differential equation /.L(x)f’(
x) +
+a*( x)f”(X)
= - 1
290 with boundary
JOSEPH C. WATKINS
AND BIRGIT
WOESSNER
conditions f(0) = 0
f(L) =o.
and
Let 7, as before, be the random time when the cell first reaches Again, we write Dynkin’s formula, ’
0 or L.
= f(x)+EQp)dr] =f(x)-E”[ps] = f(x)-
EX7.
As before, X(T) can take only two values, 0 and L. In both instances, f(X(T)) = 0 by the boundary conditions. Therefore, 0= f(X)-EXT or
Thus, T(x) solves the differential then
equation
stated above. Let g(x) = f’(x>;
j_+)g(x)+&+)g’(x)
= -1
is a first-order linear nonhomogeneous ordinary differential a’(x) is continuous and positive, then we divide u* to obtain
The integrating
factor for this differential S(y)=exp
equation
If
is
-~y$!$fdt]. [
Thus, we integrate. The constants of integration boundary conditions, and we obtain
W=@s)
equation.
can be found
dZ)dY, ic-joYu2(z;s(z)
using the
DIFFUSION
MODELS
291
FOR CHEMOTAXIS
where the constant -1
c= i ~L~(Y)~y~2(z;s(z)
dZdY)(kLSCY)dY)
.
Thus, the model and the local estimates for I_Cand a2 at the points O=&)
