.1006 UBRARY Michigan State University This is to certify that the dissertation entitled A fully explicit optimal two-stage numerical method for solving reaction-diffusion-chemotaxis systems PhD presented by Jui-Ling Yu has been accepted towards fulfillment of the requirements for the degree in Mathematics " ”I ’ ,,;.> ' w ,,‘----"“"’"r 1‘ ~ C L ”3.374411”, 1*Mv) Major Professor’s Signature 5—5-05‘ Date MSU is an Affirmative Action/Equal Opportunity Institution _ _..-.-.-.-o--n-----o-o--------o-c--'---- or— v r -- PLACE IN RETURN Box to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE A FULLY EXPLICIT OPTIMAL TWO-STAGE NUMERICAL SCHEME FOR SOLVING REACTION-DIFFUSION-CHEMOTAXIS SYSTEMS By Jui-Ling Yu A DISSERTATION Submitted to Michigan State University in partial fillfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics 2005 ABSTRACT A fully explicit optimal two—stage numerical scheme for solving reaction—diffusion-chemotaxis systems By Jui-Ling Yu Reaction-diffusion—chemotaxis systems have been proved to be fairly accurate mathematical models for many pattern formation problems in chemistry and biology. These systems are important for computer simulations of the patterns, parameter estimations as well as analysis of the biological properties. In order to solve reaction- diffusion-chemotaxis systems, efficient and reliable numerical algorithms are essential for the pattern generations. In this thesis, a general reaction-diffusion-chemotaxis system is considered and specific numerical issues are discussed. We propose a fully explicit discretization combined with a variable optimal time step strategy for solving the reaction-diffusion—chemotaxis system. Theorems about stability and convergence of the algorithm are given to show that the algorithm is highly stable and efficient. Numerical experiment results are given for one testing problem and one real experi- mental problem. To my parents, brother, sister and my grandparents iii ACKNOWLEDGMENTS I would like to express my sincere gratitude and thanks to my advisor, Professor ChiChia Chiu for her constant and patient help, encouragement, and excellent advice. I would also like to thank Professors Michael Hazier, Moxun Tang, Chang-Yi Wang, and Farhad J aberi for their time and valuable suggestions. Thanks for people in the Lansing Buddhism Association, who offer their time and space for me to learn Buddha’s teaching. Thanks for my teachers and professors Wei-Yueh Lin, Wei-Eihn Kuan, Min-Chu Lee, Charles MacCluer who give me valuable suggestions and support. Thanks for my friends, Cherry Lee, Yi-Chia Lin, Elif Seckin, Jaegwi Go, Hsiu-Chuan Wei, Kuerak Chung, Jung—Chun Liu, EunJeong Kim, Steve Gagola, Derek Lorek, Danyu Zhu, Mei-Yu Tsai, Aloka, Chin-Chieh Chiang’s family who helped me in many different ways. Thanks for my lovely cousins who supported and comforted me when grandpa passed away. Thanks for my cousin Maggie Lin. Thanks for my students, who make me become a better teacher. TABLE OF CONTENTS LIST OF FIGURES LIST OF TABLES Introduction 1 A fully explicit discretization 2 An optimal two-stage Runge-Kutta method 2.1 Stability issue ................................. 3 An optimal two-stage numerical scheme 4 Stability and Convergence 5 Numerical Experiments 5.1 1D testing problem .............................. 5.1.1 Results for 1D experiment ........................ 5.2 Simulation of Adler’s experiment ...................... 5.2.1 5.2.2 5.2.3 5.2.4 5.2.5 5.2.6 5.2.7 Chemotaxis phenomenon .......................... Experiment results ............................. Mathematical Model for pattern formation in E. coli experiment . . . . The coefficient function .......................... Intuitive explanation of pattern formation mechanism .......... An optimal two-stage numerical scheme ................. Results of numerical experiment ...................... 6 Conclusion BIBLIOGRAPHY vi viii 13 15 20 27 31 37 38 39 41 46 50 51 59 63 65 1.1 2.1 5.1 5.2 5.3 5.4 5.5 5.6 5.7 LIST OF FIGURES Method of Line ................................ The solution of the MinMax problem ..................... Error between numerical simulation and analytic solution. ........ Photographs showing bands of E. coli in a capillary tube. In all the exper- iments reported here, capillary tubes were filled with a liquid medium, inoculated at one end with 200000 to 2000000 bacteria, and then closed at the ends with plug of agar and clay, all according to procedure de- scribed in full elsewhere. The tube were incubated horizontally at 37C. The origin, which is turbid because of the bacteria that have not moved out, is visible at the left, then the second band of bacteria, then the first band. Plugs at the ends are not shown. The concentration of galactose was 0.00025 mole per liter. .................. Utilization of oxygen by bands of bacteria in 0.00025 M galactose. Oxygen was measured polargraphically by inserting an oxygen needle electrode into the capliiary tube in 4—mm steps when the first and second bands were visible where shown by arrows .................... Utilization of 0.00025 M galactose. The tube contained 014- galactose (0.00025 mole/liter and 1000000 counts per minute per milliliter). At 4 or 8 hours, when the first and second bands were visible where shown by arrows, the tube was factionated into ten compartments, each 8 mm long. the contents of each compartment were chromatographed on paper, with n—butanol, acetic acid, and water (12:3:5) as the solvent, and the radioactivity in the galactose region of the chromatogram was measured in a paper-strip counter. ................... Utilization of oxygen by bands of bacteria in 0.000005 M galactose. Mea- surements were made as in Figure 5.3 and Figure 5.4. Utilization of 0.000005 M galactose. The tube contained Cl4—galactose (0.000005 molde / liter and 200000 counts per minute per milliliter). At 4.75 hours 9 18 38 41 43 when the first and second bands were visible where shown by the arrows. 45 The relation between the amount of galactose and the distance of the bands of bacteria .............................. This set of pictures simulates the movement of bands of bacteria when the concentration of galactose is higher than oxygen ............. 59 60 5.8 The picture shows the movement of bacteria when the concentration of galactose is below the oxygen. ...................... 61 5.9 Relation of the amount of oxygen and the distance of the bands of bacteria. 62 vii LIST OF TABLES 5.1 Number of iterations .............................. 37 Introduction One of the important branches in mathematical biology is the study of pattern for- mation. In 1952, Alan Turing proposed the reaction-diffusion model as the chemical mechanism for biological patterns [25]. The idea was that the interaction between the passive diffusion and chemical reaction will envolve into a spatiotemporal pattern. Moreover, the motion of the biological cells is often influenced by the environ- ment. This influence refers to chemotaxis, which is often known as directed random movement. The presence of chemotaxis provides a general tendency for bacteria to move toward the chemotaxattractant [17]. The study of chemotaxis phenomenon among bacteria began in the nineteenth century. In 1884, Wilhelm Pfeffer, Engelmann, and their associates discovered that bacteria move preferentially toward a concentration of oxygen, minerals and organic nutrients. These researches demonstrated chemotaxis by observing whether bacteria in a suspension gathered near or away a gas bubble which contains a certain amount of oxygen [4]. In 1893, Beijernck discovered that bacteria move toward oxygen “ macroscopically. He concluded that the bacteria search for a certain optimum concentration of oxygen.” Later, Sherris, Preston, and Showmith and Baracchini and Sherris comfirmed the results. A significant number of researches in molecular biology were done by Julius Adler and his co—workers around 1966. He tested the chemotaxis behavior of bacteria by using Escherichia coli and Salmonella typhimurium bacteria. His strategy of studying the chemotaxis phenomenon by using the simplest organism has provided the proof for and inspired the study in the mechanism of synthesis and control of biology macromolecules. In this thesis, we studied one of Adler’s papers on the chemotaxis behavior of the Bach bacteria because “ the vast knowledge of its biochemistry and genetics could be brought to bear on the problem. [18]” The mobility of the E. coli and many other bacteria depends on swimming and tumbling. Swimming is a near straight run. Tumbling is a very rapid and jerky movement during which the cells are reoriented to a new direction. The movement of swimming (lasts for seconds) is interrupted by tumbling (lasts for tenths of seconds). Cells respond to the temporal and spatial changes in the environment by modulating the probability of tumbling. When cells encounter a favorable environment, they reduce the frequency of tumbling so that the swimming time in a favorable direction is increased. Therefore, a net migration along the desired direction is formed. That is chemotaxis [8]. The motion of swimming and tumbling is controlled by bacteria flagella which consists of a long filament. It is a left handed helix. The flagella turns either clockwise or counterclockwise as a propeller. The clockwise motion producing a pushing force results in the motion of swimming. The counterclockwise motion resulting in deformation of the filament produces tumbling [8]. Bacteria sense the change in the local concentration of chemicals through membrane receptors binding chemical’s moleculaes. The cells estimate the change in the gradient of chemicals by measuring the change in the receptor’s occupation [16]. It results in a shift in the clockwise and counterclockwise motion in the flagella motor. The addition of an attractant results in the supression of tumbling. The supression time may last for a few seconds or minutes. It depends on the strength of the chemotaxis [8]. The general form of the reaction-diffusion—chemotaxis system is % = DuAu _ V(ux(a)Va) + f(z, y. t) a a“ = V - DaVa + 90w) u(:r, y, t), a(:r,y,t) represent the concentation of cells and the concentation of chemoattractant substrate respectively. DuAu is the diffusion term or random walk, and which mainly characterize the movement of swimming. Du > 0 is the constant diffusivity describing the speed of bacteria. V(ux(a)Va) describes the biased random movement which is induced by external chemicals. It has been proven that the nonlinearity in the reaction term is a necessary condition for the pattern to be formed. No general analytical solutions are available at the present time. Thus, a numerical simulation provides a comprehensive understanding of the reaction- diffusion-chemotaxis systems and the mechanism of the underlying biology [17]. For the past years, numerical methods for reaction-diffusion systems have been well-studied: for example, the ADI type of scheme, the finite element method, monotone iterative methods, and the multigrid method. However, the development of numerical methods for solving reaction-diffusion—chemotaxis systems have not been thoroughly studied due to the additional nonlinear chemotaxis term [26]. Furthermore, most of the time, numerical methods are too complicated to be used 3 in practice. The idea of this thesis is to construct a numerical method, which is not only reliable but also easy to implement so that it can be applied to complex systems. Many numerical methods have been studied to solve the PDEs: for example, the fully implicit method, the ADI method, and the explicit method [3, 23]. The explicit method is especially easy to implement and does not have to invert the matrix resulting from the implicit method. One of the approaches to get the explicit scheme is to use the method of line, which discretizes the spatial and temporal variables separably, changes PDEs to an ODE system, and solves the ODE system at each time level [1]. In this thesis, the standard finite difference and the upwind scheme are proposed for spatial discretization, and the spatial matrix obtained from it are sparse matrices. They are less expensive to integrate than the full matrix during the numerical simulation [11]. The convergence rate of the numerical solution strongly depends on the distrib- ution of the eigenvalues of the spatial discretization. If the matrix is symmetic and positive definite, many efficient iterative methods can be achieved to obtain fast convergence. The spatial matrix we obtained from the finite difference scheme is not symmetric due to the chemotaxis term and the boundary condition. Therefore, it is necessary to develop an efficient numerical scheme to treat the nonsymmetric case. The development of numerical methods for solving the ODE system are well— studied: for example, Euler’s method and the n‘h—stage Runge-Kutta method. The explicit Runge-Kutta method has advantages of low storage, high accuracy, flexible stability region and easy implimentation [11]. However, since the method is explicit, the time step is often restrictive. It has become an unfavorable feature of the explicit scheme. There are many research papers addressing the issue of improving the convergence rate for Runge—kutta method with different stages [11, 12, 13]. Since the spatial matrix of the ODE system is time dependent in the chemotaxis system, the stability region will vary by time. Since the size of the matrix is very large and the matrix is not symmetric, it is not possible to find the eigenvalues analytically, and using the existing numerical methods to estimate the distribution of eigenvalues are too expensive. Therefore, the time step has to be forced to be very small to satisfy all the stability conditions for the ”time dependent” ODE system. In this thesis, we develop an efficient numerical method. The anar- lytic form of explicit time can be easily calculated and the time step is locally optimal. We study the properties of a particular Runge—Kutta method with n=2 for choosing the optimal time step. Since for nth-stage Runge-Kutta method, the stability function is a nth degree polynomial function, the analysis for the MinMax problem becomes complicated and the analytic form of the optimal time step does not seem possible. Furthermore, introducing more stages does not guarantee the efficiency of the CPU time [12]. The generation of patterns depends on the distribution of eigenvalues of the ODE system. The ODE system contains eigenvalues with positive and negative real parts. We have to take these two parts into consideration during the numerical simulation. The stability issue has to be addressed since we want to make sure that patterns which are generated from the computation simulation are derived from the original PDEs, not due to numerical instability. Ian. We propose an optimal two-stage method to solve the ODE system with time dependent spatial matrix. The method keeps the advantage of the two-stage Runge-Kutta method for its higher accuracy and easy implementation. In addition, the time steps are enlarged, also we can arrive at the explicit formulation for the time step. The numerical experiments show that the optimal two-stage scheme is accurate and converges faster than the standard two-stage Runge-Kutta method for the time evolution problem. In Chapter 1, the fully explicit discretization for the spatial variables is discussed. The method of line and the upwind scheme for the chemotaxis term are the main topics of this chapter. The method we used to find the optimal time step for the standard two-stage Runge-Kutta method will be addressed in Chapter 2. In this chapter, we mainly discuss the method to solve a MinMax problem. In Chapter 3, a fully explicit optimal two-stage scheme is proposed to solve the ODE system. The consistency of the numerical method is investigated. Stability and Convergency are addressed in Chapter 4. Finally, the results of numerical experiment on the comparsion between the standard two-stage Runge—Kutta method and our optimal two-stage scheme will be shown. A real problem will be tested against our numerical method. CHAPTER 1 A fully explicit discretization Many interesting patterns happen as bacteria consume nutrients, proliferate, perform random walks and the biased movement towards high concentration of chemicals [4, 10, 16, 21]. The time evolution of chemicals is described by reaction—diffusion- chemotaxis equations [20]. In these models, the microorganisms are simulating through 2D density, and a reaction-diffusion—chemotaxis equation of this density characterized their time evolution. The equation is coupled with other reaction- diffusion equations, since the chemotaxis substrate diffuses and produces by itself [16, 20]. The general reaction—diffusion—chemotaxis systems have the following form: % = DuA'u. — V(UX(G)Va) + f(xa y, t) 50. E _ v.D,Va+g(a,U)- u(:r, y, t), a(x,y,t) represent the concentation of cells and the concentation of chemoattractant substrate respectively. DuAu is the diffusion term or random walk, and which mainly characterize the movement of swimming. Du > 0 is the constant diffusivity describing the speed of bacteria. V(ux(a)Va) describes the biased random movement which is induced by external chemicals. The movement is biased either along the gradient direction or in the opposite direction of the gradient [16]. x(a) 2 0 is the chemotaxis parameter function describing the speed of chemotaxis, f(:1:, y, t) is the reaction term for cells, which indicate the growth rate of cells. For simplicity, we use a 1D reaction-diffusion-chemotaxis equation to demonstrate the analysis of our numerical scheme. We also assume boundaries are periodic with period 3, and the initial value is known, Bu 6—1; = uuza: _ (UX(G)az)x + f($,t),.’L' 6 [0,6], (1.1) u(z,0) = u0(:r),x E [0, Z], u(0, t) = u(€, t). The equation (1.1) is generally nonlinear so the analytical solution is difficult to solve. Therefore, numerical simulation becomes crucial to investigating the mecha- nism of the underlying biology. To solve this equation numerically, a very common method is to apply the method of lines [15] which discretizes spatial and temporal variables separately, then turns the PDE into an ODE system and solves the ODE system at each time level. Take a uniform mesh of f2 = [0, K] with mesh size h = Z/N. Let 1:,- E Q be a grid point, then x,- + 1 = 2:,- + h. A time step of size t" + 1 — tn will be denoted as At, and we adopt the standard notation a? z u(:r:,~,tn), etc. We use the standard second-order centerd difference to approximate the Laplace operator, and 8 Figure 1.1: Method of Line the first-order upwind scheme for the gradient operator. The assumption of using an upwind scheme [19] is crucial in the deveopment of our numerical method. It ensures the stability of the numerical scheme. For simplicity, let us denote x(a)Va by q(t). Then, Bu 5 = dunm — (q(t)u)z + f(a:, t) = Muir... —2ur+ur_1 h2 )_ (qzu + quit)? + fin, where (gum)? is approximated by the first-order upwind scheme, is. q”+ q” n (qua? z ‘ 2,] "van.- (q. > 0) If— I‘ + ——q1 2hlqzl‘711ui, (q? < 0) (1.2) (1.3) (1.4) where ‘77-f-a:ui = uil+1 — u? If q? > 0, (1.2) becomes 0a 2a" +u" q" 5,; = duth‘z—) -( an)” - flu: - at: 1) + f." du n —2du — h—Z—ui+1+(T—(qz)"—q—Wll +(% )11+fin(15) +h The periodic boundary condition can be approximated by UN—i = “—1 UN+1 = U1 Let u" = [11?] and f” = [fin] be a vector obtained by the usual ordering. Then (1.5) can be represented by the matrix form as follows: = M(t")un + f". (1.6) d— Letwi=—_ig‘u- (9x)? ‘gfi-Z 0,1,2,---N, then M(tn)= 10 \ E O «as are + s|f= _/ / 3:35“ O 3345‘ + SE: am} a (N+1) x (N+1) matrix. Observe that the spatial matrix is time dependent. By the method of line, we have turned the PDE into an ODE system. Many numerical methods have been studied to approximate the temporal derivative in the ODE system. For example, Euler’s method, and the 72‘” stage Runge-Kutta method [15]. Among all of these numerical schemes, the explicit Runge-Kutta method has been widely applied. It has the advantages of high accuracy, flexible stability region, low storage, and ease of implementation [11]. However, step size has been limited due to the nature of the explicit method. Due to this fact, we are forced to choose a smaller step size than is needed to guarantee the convergence of the numerical scheme. Oftentimes, the restricted time step increase the CPU time which becomes a unfavorable feature of the explicit scheme. Many researchers have been working on improving the preformance of explicit Runge—Kutta method by releasing the time 11 step constraint. Our idea is to adjust the appropriate time step so that the spectral radius of the spatial matrix in the ODE system is minimized. Therefore, optimal convergence can be achieved [11, 13]. This leads to the solution of the MinMax problems [11, 13]. Detailed analysis of the MinMax problems will be addressed on Chapter Two. CHAPTER 2 An optimal two-stage Runge-Kutta method For the nth stage explicit RK method, the stability function 12(2) is a function with degree 11. When n 2 3, the formulation of R(z) becomes complicated and the analytic form of the optimal time step does not seem possible. Furthermore, it is too expensive to update the time step by using the modified projected Lagrangian Algorithm and introducing more stages of the optimal RK method does not guarantee the efficiency of the CPU time [12]. Therefore, We focus on applying the optimal two-stage RK method to solve the chemotaxis system under the stability consideration. The main purpose of this section is to obtain the optimal time step for the standard two-stage Runge-Kutta integration. First of all, let us assume the spatial discretizar tion matrix M and the reaction vector in the ODE system is time independent. du 13 If the we set the reaction vector as zero, the system becomes du_ Let us use the standard two-stage Runge—Kutta integration to solve (2.2), we get 11(0) = u" + AtMu” (2.3a) an“ = a" + 92i(Mu" + Mum). (2.3b) The parameter At usually stands for the time step. When we treat (2.3) as a stande iterative method, it is actually a one-step method: A u”+1 a" + %Mu" + 7tM(un + AtMu”) At2 2 n ll Let P denote the amplification factor, i.e, Ar2 2 then (2.4) can be written as 11"“ = P(At)u" P"+1(At)u°. (2.6) The generation of patterns depends on the distribution of eigenvalues of the ODE system. The system contains eigenvalues with positive and negative real parts. 14 We have to take care of these two parts during the numerical simulation processes. Therefore, the stability issue has to be addressed. 2.1 Stability issue Let us consider the nonhomogeneous ODE system, du n n where M is a constant matrix, f” is a constant vector. Let v0 denote the perturbation of the initial condition no. Then the solution vector v” + 1 after nth iterations and the exact solution a" + 1 satisfy the following: n+1 n. At2 n u = P(At)u +(1+At+—2—)f (2.8) n+1 n Atz n v = P(At)v + (1 + At + T)f . (2.9) Let 6" + 1 = v" + 1 — an + 1 denotes the perturbation error after n + 1 iterations, we have 5"“ = p(At)e" = p"+1(At)5°. (2.10) Horn equation (2.10), we know that the error e" + 1 will be bounded if the norm of the amplification factor is less than one. Nevertheless, the iterative procedure will converge, if and only if, the spectral radius of the amplification factor is less than a unity (Issacson and Keller 1966). Moreover, the convergence speed of the iterative 15 procedure depends on the size of the spectral radius; if the size of the spectral radius is small, the iterative procedure will converge faster. This gives rise to formation of the MinMax problem. The characteristic polynomial of the amplification factor is er A(At,/.t) = 1 + Atp + 2 #2, where a is an eigenvalue of M. Note that the characteristic polynomial A(At, a) is a function of time step and eigenvalues. Definition 1 Let a = a + bi. Define f(a, b, At) = |A(At,p)[2. Theorem 2 If b is bounded, then f(a, b, At) = f3(At) = (1 + aAt + %a"’At"’)2 + 0(At3). Proof. 2 f(a, b, At) = |1+ At(a + bz‘) + ATtnz + bi)2|2 - Atz 2 2 -2 |1+At(a+bt)+7(a —b +2ab2)| 2 [(1 + aAt + ATthzZ — b2)]2 + (Atb + Atzab)2 4 = (1+ aAt)2 + Arm;2 — b2)(1+ aAt) + ATtm? —— (22)2 (PM + 0(At3) + 4 (1 + aAt)2 + a2At2(1 + aAt) — PM + ATta“ (2.11) + b2At2 + 0(At3) + 0(At4) = (1 + aAt + gazAt2)2 + 0(At3) (2.12) Since the truncation error of the standard two-stage Runge—Kutta method is 0(At3), the higher order terms appeared in can be ignored. Moreover, the min At maxa f (a, b, At) can be simplified as min At maxa fa(At) since we are only interested in the case where % S fa,(At) < 1 when a < 0. It is worth noting that the optimal two-stage scheme matches well with the theoretical analysis in the following sense: eaAt is numerically represented by fa(At) = 1 + aAt + §a2At2. If a > 0, the accuracy of the numerical solution is the only criteria for us to choose the time step. We have to choose At small enough so that fa(At) will accurately integrate eaT. The stability issue is not important here, since eaT is growing and aAt will not lie in the stability region. If a < 0, we have to consider both accuracy and stability of the numerical scheme. However, they often imply each other. When a < 0, fa(At) corresponds to a decay mode. We have to choose the right time step so that fa(At) is less than a unity. In other words, the choice of the time step has to restrict aAt to stay in the stabiliy region so that fa(At) is a correct integratration of eaT[15]. Lemma 3 For o < 0, fa has the following properties: (a) fa(0) =1 and rot—72) = 1. (b) fé(0) = a. (c) fa(At) 2 0, At > o, and fa(At) < 1,At e (o, —-2/a). (d) mom 6 R fa(At) occurs at At = —1/a. (e) minAt e R fa(At) = 1/2- Proof. The proof is straightforward using the properties of the second—degree poly- nomials. Cl Theorem 3 shows that the optimal time step can be determined by the largest and smallest real part of the eigenvalues of the ODE system. Theorem 4 Ifa < 0 and -amax < a < ‘a'min’ where “min > 0 and amax > 0 are constant then min At > 0 maxa E l—amax ‘am' ] Ifa(At)|2 is obtained at Atopt = i in 2 amin amax' Proof. f;(At) = a + a2At = 0. Let f;(At) = 0, get At = —%. Moreover, f(—$) = g > 0. ram 1 NIH llaMax Alopt -1/a 1/aMin Figure 2.1: The solution of the MinMax problem. From Lemma (3) and the graph, it is clear that the optimal time step will be obtained if f—amin(At0pt) = f—amax(At0pt)- 18 i.e. 1 2 2 1 2 2 1 — aminAtopt + 5amnAtopt = 1 — amaxAtopt + -2—amxAt opt 5(amax - airin) Atom = amax _ amin- which implies 19 CHAPTER 3 An optimal two-stage numerical scheme In this section, we will extend ideas from chapter 2 to the nonhomogeneous ODE systems with time dependent spatial matrix, i.e, du d? = (t)u + f(t). (3.1) An optimal two—stage scheme will be constructed. The consistency of the numerical scheme will be discussed in this section. Let us consider the following two-stage scheme for solving (3.1) ”(1) = u" + alAt[M(t") + M(t"+1)]u" (3.2a) on“ = u" + agAt[M(t") + M(t"+1)]u(1). (3.2b) There are two parameters a1 and a2. We are looking for a} and a2 so that the two-stage scheme is consistent in time and space, i.e, having the truncation error 20 being 0(At3). Meanwhile, the time step will be locally optimal according to the analysis in Chapter 2. Treating (3.2) as a standard iterative method, the two-stage scheme becomes un+1 ll u" + agAt(M(t”) + M (t"+1))[u" + alAt(M(t") + M (t"+1))u"] [I + a2At(M(t") + M (t"+1)) + a1a2At2(M(t") + M(t"+1))2]u"_ Since Mtt" + 1) = Me") + AtM'rtn) + om). u"+1 [I + azAt(M(t") + M (t") + MAW”) + 0W2)” + a1a2At2(M(t") + M(t”) + AtM’(t”) + 0(At2))2]u" [I + 2a2AtM(t”) + o2At2M'(t") + 40102At2M(tn)2 + 0(At3)]u". (3.3) Moreover, Taylor expansion of u at time t” is o(t"+1) = u'(t”)+u(t")At+ WT )At2+0(At3) = u(t") + AtM(t")u(t”) + -A2—t(M’(t")u(t”) + M(t")u’(t"))+0(At3) = u(t") + AtM(t")u(t") + A—2t2-(M’(t")u(t") + M(t")2u (t”))+0(At3) = [I+AtM(t")+A thM’n") + A—t2M(t")2]u(t") + 0(At3). (3.4) 21 Matching coefficients from (3.3),(3.4) get a _1 1—4) and 1 2. C12: Therefore, we can set up the optimal two-stage scheme as the following: um = u” + iAt[M(t”) + M(t"+1)]u" (3.5a) on+1 = u" + éAt[M(t") + M(t"+1)]u(1). (3.5b) Note that the coefficient of the Euler’s predictor is 1, and the coeflicient of 'ltapezoid corrector is % [15]. (3.5) is not a standard two-stage Runge—Kutta method. The con- sistence of the numerical scheme will be derived in Theorem 5, which shows that the difference between the differential equations and the two-stage scheme we presented here can be controlled. Because the analysis for the 2D case is similar to the case for 1D, for simplicity, the 1D model equation is presented in the proof. Theorem 5 (Convergence) Assume that ][ Mn+Mn + 1 [IS a for alln, let |] en ||=|] on — on") n, and (n + 1)At g T. If || 60 u: 0 then u e" “g 0(h) + 0(At2). Proof. The general reaction—diffusion-chemotaxis eauation is an a = duuzx _ (q(iL‘)U.);,-, + f(IE, t)’ (3'6) 22 Since the upwind scheme is of order one in space, the finite difference discretization gives Bu _ uf+1—2u?+u[’_1 git—Gill , 5i _ d“( h2 l h 1“ q"+l<1”| qt-lqfl n _ 1__1. n ui_;__1_ n i . Oh 2h v—zu 2h V+xu +f1. + ( ) du _2du = h_2 ui+l+(— [7,2 qt _qi— _qi + qt 11 qt" n 61‘“ n — ,1 2,] "'1. — —-2—q"",,' 'V.u.> +h3u._.-h1+f+0() wheret=0,1,2,---,N. Therefore, an n n n 5Y=M(t )u +f +O(h). (3.7) Applying the two-stage scheme to (3.1), we get um 11" + iAt[(M(t") + M(t"+1))u" + f(t”) + f(t"+1)] (3.8) on“ = a" +éAt[(M(t")+M(t"+1))u(1)+f(t")+ f(t”+1)] (3.9) Viewing two-stage scheme as a standard iterative method, we get n+1 n n+1 11"“ = [I+At(—M +M +E(———M +M )Zun 2 2 2 + _[(fn+fn+1)+_ 4t(Mn+Mn+1)(fn+fn+1)] The Taylor expansion of u at t" is u_(t"+1) — u(t")+u’(t")At+u "Us —)At2+0(A At3) 23 = u(t") + At(M(t")u(t") + f(t") + 0(h)) + ATttM'unuu") + Means") + Maura") + f(t")’ + 0(h)) + 0(At3) Since M(t"') = M(t" + 1) — M’(tn)At + 0(At2), U( tTH-l) = u(t") + %M(t")u(tn) + %(M(t"+1) — M’(t")At 0(At2))u(t") + At f" + At0(h) + ATtZ(M’(t")u(t") M (t")2u(t") + M (t") f" + (f")’ + 0(h)) + 0(At3) 2 u(t") + %(M(t”) + M(t"+1))u(t") — 92LM(tn)'u(tn) 0(At3) + At f" + At0(h) + -A2i2(M’(t")u(t") + M (t”)2u(t")) 2 2 2 A7tM(t")f” + AT’W + AT’OUL) + 0(At3) At2 2 2 2 _M(t")f" + At f" + At0(h) + —2—0(h) + —; (f")' 2 2 12 ac") + 923(Mlt") + M(t”+1))u(t") + ATHM“ ) u(t") + %(M(t") + M (t"+1))u(t") + M (t")2u(t") + 0(At3) 2 {@222 _ _M_§‘”_)At + 0(At2))2u(tn) + 0(At3) 2 2 ATtZM(t")f" + Atf" + At0(h) + -A-2t—0(h) + ATtUn), (I + 923w" + Mn“) + é§(M(t") + M(t”+’)2)u(t”) + 0(At3) é2fiM(tn)fn + Atf" + é;(fn)l + At0(h) + ATROUL) Assuming there is no rounding error, i.e, 24 f" = f(tn),M(tn) = M”, and en+ 1 = u(tn+ 1) — un+ 1, then n+1 2 e = [I + 925(M” + W“) + —A8t (M" + M"+1)2]e" + 0(At3) At2 At2 At2 + 7M"? + Atf" + 7f; + At0(h) + —2—O(h) At At _ 7Un+ fn+1+ T(jwn + Mn+1)(fn +fn+l)] 2 = [I + %(M" + M"+1) + A%(M" + Mn+1)2]e” + 0(At3) + At0(h) At2 At2 At + 7001) + Atf" + TU"), — ”2—(fn + 1m“) At2 At2 + TMnfn—?(Mn+Mn+1)(fn+fn+l) At n n+1 At2 n n+1 2 n 3 : [1+7(M +M )+?(M +M )]e +0(At)+At0(h) At 2 + Atf" + A715(f")’ — £22m — —2-(f" + (f")’At + 0(At2)) At2 At2 + —2—M"f" — ?(2M" + O(At))(2f" + O(At)) At n n+1 Atz n n+1 2 n 3 = [I+?(M +M )+?(M +M )]e +0(At )+At0(h) 2 = [I + %(M" + M"+1)+ fist—(Mn + Mn+1)2]"+1e0 + (n + 1)At(0(h) + 0(At2)). Since (n + 1)At S T and co = 0, it implies ||e"|| g 0(h) + 0(At2). This completes the proof. I] n n + 1 Let Jn = M t + M t , then the two-stage scheme for the homogeneous problem can be written as 25 u" + 1 = (I + AtJn + éAtnght". It follows from the analysis in chapter 2 that the optimal time step will be obtained if we set Atopt = W. Let {aili = 1 be the real parts of the eigenvalues of Jn, then “min = mz’nai < 0 {at} 6'mass = "ma,- > 0 {02'}- Observe that the matrix I + AtJn + éng depends on time. The stability region varies as the numerical scheme is marching on time. In practice, the optimal time step will be updated according to the distribution of the eigenvalues of the matrix I +AtJn + 9%ng at each time level. We call the two-stage scheme an optimal scheme in this sense. 26 CHAPTER 4 Stability and Convergence In this section, we will investigate the stability and convergence of the optimal two-stage scheme. Theorem 6 (Convergence) Assume that [I Mk + Mk +1 [[3 a for all k + 1 5 N, (tN S T, tN+1 > T) and tn = nAt g T. If ]|u(t0) — 110]] = 0, then I] u(tn) — n [[5 C(T)(0(h) + 0(At2)). Proof. Let ”en” = ||u(tn) — an”. From the proof in Theorem 3, we get n+1 At Atz 2 n 3 lie 1 s II1+7+YHue “+0(At) + At0(h) 2 3 (1+ %a + égt—azfllenll + 0(At3) + At0(h) 2, g (1 + %a + AYE?)2 n e"‘1 n At At2 2 3 3 + (1 + 70: + ?a )(O(At )+ At0(h)) + 0(At ) + At0(h) 2 g (1 + %a + est—ozrfl u e0 n 27 At At2 2 n At At2 2 ,,_, + [(1+—2—a+?a) +(1+—2—a+-——8—-a) + ---+ 1](0(At3) + At0(h)). Since ”.30” = o, and At At2 2 a At 2 _ __ = _ __ A 1+2a+8a 1+(2+8cx)t 3 1+ fiAt, therefore, A 2 . (1+ 95th: + 73:02)”+1 _<_ em. Moreover, since 2 a a (I — —At > —, 2 + 8 _ 2 it implies that 1 2 a 0:2 S _' 5 + §At CY We conclude that n (1+ S‘ia + Muzzy“ — 1 || 6 +1 H = 292% +8-ég3-a2 (0(At3) + At0(h)). g EefiT(O(At2) +0(h)) C! < C(T)(0(At2) +0(h)). U Lax’s equivalence theorem says that if we approximate a linear initial-value problem by a linear finite-difference method that satisfies the consistency condition, 28 then stability is necessary and sufficient condition for convergence [23]. Since the reaction-diffusion-chemotaxis equation is nonlinear, it does not satisfy Lax’s condition. Therefore, we have to investigate the stability issue. Theorem 7 (Stability) Assume that H Mn + Mn +1 “_<_ a for all n + 1 S N, let [en |= marl | “l —— u(tl) [,0 $1 3 n for all n such that nAt _<_ T. If] e0 [S 5 then I 6" [_<_ Cl(T)€ + C (T)(O(At3) + AtO(h)), where C = C(T) is a constant which depends on T, but not At and Ax. Proof. It follows from the proof of convergence, we have ,, At ll 6 +1” S H 1+ -2—(M,, + Mn+1) At2 + —8—(M.. + Mn+1)2 II H e” II +0(At3) + NOW 3 (1+ 923a + éga’) ll 6”“1 II + (1 + 921a + %fia2)(O(At3) + AtO(h)) + 0(At3) + AtO(h) 3 (1+ %a + %Q2)n+15 + [(1+ %a + 953022)" + (1+ %a + %fia2)"“1+~-+1](0(At3)+ AtO(h)). Since (1+%a+%§3a2)n+(1+%§a+-‘}2fia2)"’1+---+1 29 therefore, (1 + %a + —Z%3a2)”+15 + C(T)(0(At3) + AtO(h)) (1 + fiAt)"+le + C(T)(0(At3) + AtO(h)) ll 6”“ II |/\ |/\ o16o + C(T)(0(At3) + AtO(h)) |/\ l/\ e37}: + C(T)(0(At3) + AtO(h)) l/\ efiTe + C(T)(0(At3) + AtO(h)) S 01(T)5 + C(T)(0(At3) + AtO(h)), where Cl (T) = efiT and C(T) = 2/ozefiT. El 30 CHAPTER 5 Numerical Experiments In this chapter, we use a 1D testing equation to test the efficiency of the two-stage numerical scheme. The results of a comparison between the numerical and analytical solutions will be given. Furthermore, we will show that the optimal two-stage scheme converges with the analytical solution for the time evolution problem while the standard two-stage Runge-Kutta method converges only when the time step is small enough. Moreover, the optimal two-stage scheme converges much faster than the standard two-stage Runge—Kutta method in the region of convergence. Later in this chapter, we will test the numerical scheme against a realistic problem. 5.1 1D testing problem Let us consider the following testing equation du where f = (ta: — 1)e—tcos(:r — t) + te_tsin(a: — t). 31 Assume :t: G [0, 2n],t > 0, and the boundary is periodic. The analytical solution of (5.1) is given by t u(a:, t) = e‘ sin(a: — t). We solve the equation by the two-stage scheme. We use the first-order upwind method to approximate the chemotaxis term, and the standard second—order centerd difference to approximate the diffusion term. Since t2: is positive for t > 0, the chemotaxis is approximated by the backward difference operator. This leads to the following setting: Let u2-(t) z u(:r2-,t), then du “4+1 '— 2114' + ut—l ui — rut—1 Et- : h2 —tu2 —t£L’i( h )+f($i,t) 1 -2 tn 1 txi == mum + (723' — t -- ‘h—lui + ('2'; + ”Ii—)ui-l + “$13 t)’ where i = 0, 1, 2, ...N. The boundary condition is u_1= ’LtN_1. Therefore, we can set up the ODE system as the following : du_ __Mn 11 d, (t)U+f. where M (tn) = 32 n 1 tux] 1 _ fl \ f TO B? "l' h it? h 1 tna: 1 '22 + —hQ T? ‘52 1 tnx n 1 + 7' i? 2 h. 1 1 t"a:N_1 l a a + h it i where T? = TI; —tn— 1%,,- =0,1,2,~- ,N. _ n n + 1 The matrix M n = M“ ) + 340 ) can be evaluated by changing the time level from t" to t" +1. Since the matrix size is usually very large and entries of the matrix vary according to different time levels, it is not possible to find all the eigenvalues in analytic form. The existing numerical algorithms, e.g. the power method, are too expensive to estimate the eigenvalues. Moreover, we only need the extreme eigenvalues to calculate the adapted time step. The Gerschgorin Theorem from linear algebra will be used for this purpose. Theorem 8 (Gerschgorin) The union of all discs K2- 2: {u E l/t — a22- [g 22 = 1 k 72 2. | aik I} contains all eigenvalues of the n x n matrix A = [aikl- 33 Proof. See the reference [23]. CI The estimation of the eigenvalues for the first row and the rest of the rows will be discussed separately, since they have different structures which are due to the effects of boundary conditions. For the ith row, where i = 1,2,. . . N, the diagonal entries are _ 7',-"’+'r1"+1 _ —2 (t"+t"+1) x,(t"+t”+1) — 2 — h? 2 2h ° an The absolute sum of the off diagonal entries for each row is i [0" l‘ 2 + :1:,-(t" +t"+1) . . . 2’ _ h2 2h ° 2:00752 Gerschgorin’s Theorem implies all the eigenvalues of M are contained in the union of —2 t" + t"+1 at" + t"+1 2 a:,(t" +tnt1) /\ — —- — —— — . I (h? ( 2 l ( 2h. - h2 + 2h The smallest and largest real part of the eigenvalues satisfy tn + tn+l ( tn + tn+1 $2,011 + tn+l) 2 2 ) h —4 )SASfi—( For the extreme case, _—4 t"+t"+1 :c,(t"+t"+1) ‘amin-n‘l—t—l‘ h n n + 1 Since —amax = —(t—+—ttz—) < 0, it shows that the largest eigenvalues are on the left hand side of the complex plane. Therefore, we only have to consider the case for -ama:r < 0. By theorem 4 and 1 1 AtrowO : 2 tn + tn+1 > Atrowl z 2 tn + tn+1 (tn + tn+1) > 715 + ( 2 ’ E + ( 2 ) + h h At 2 = 1 > Tow 2 t" + tn+1 (tn + tn+1) }? + (—2—) + 2h h At = 1 > rowN 2 tn + tn-l-l (tn + tn+1) ’ E5 + (—2—_) + W h the time step is 2 1 Atom = amin + amax : 2 t" +13"+1 (t” + tn“) ' (52) it? ( 2 ) + W h For the first row of the matrix M, i = 0, a _ 2 t” +t"+1) 00 — hz 2 i 1 tn + tn+1 1 tn +tn+1 Ejillaod- [ZW—l-T‘l’ I iii—(Tll The location of eigenvalues for this row is 2 tn+tn+1 1 tn+tn+1 1 tn+tn+1 IA+E§+(—2—)ISE+(——2—)+l§—(Tll- 35 In this case, —3 1 tn+tn+1 1 1 tn+tn+1 ——t" tn+1————-——— Ihg < 2h )1- _ h2+|h2 2h )I The lower bound is _3 1 tn+tn+1 _ min=__ t" tn“ _ __ ..____, a h2(+)|h2(2h)l The upper bound is 1 1 1t"+t"+1)l m" h2 lh2 ( 2h When —amax Z 0, we set -amax = 0. Therefore, the time step is separated into two cases as follows: 1 if?“ < tn+tn+l 4 + (t“ + t"+1)(4h2 + h) h2 — 2 2 2+tn+tn+1 lffi> 2 h2 2 Since time step in (5.3) is greater than time step in (5.2), the optimal time step is 1 Ato = . 5.4 Pt 2 tn + tn+1 + (tn +tn+1) ( ) PM 2 )fl 2 36 5.1.1 Results for 1D experiment Method / Time 1 2 3 4 5 Optimal two stage 284 637 1060 1553 2116 Standard Runge-Kutta At = 0.002 500 1000 1500 2000 2500 Standard Runge—Kutta At = 0.004 250 Div Div Div Div Table 5.1: Number of iterations. This table is the number of iterations with different methods at different time. The table says that the standard two-stage Runge-Kutta method will converge when the time step is small enough. However, the optimal two-stage method converges faster than the Runge—Kutta method when both methods converge. When the time step for the Runge-Kutta method is increased, it will diverge. 37 Optimal two-stage vs. Runge-Kutta two-stage method § 3.5 ~ - 4} 3 - - Two stage -1 an Ringo-Knits I A n. a g 2'5 h at '4 c 0 g a. 2 v _ .E * ‘5 it t 1.5 - , ~ Lu 1 * 1 ~ Q -1 Q 0. fl» 0 3. * 05F: \ _ P i i 0 1 3 4 5 6 7 Time Figure 5.1: Error between numerical simulation and analytic solution. This picture shows the relative error between numerical simulation and analytic solution. While the two—stage Runge—Kutta method blows up in a short time, the optimal two-stage method is still very accurate. It shows that the optimal two-stage scheme not only is accurate but also converges faster than the two-stage Runge-Kutta method. 5.2 Simulation of Adler’s experiment In this section, we will test the optimal two-stage scheme with a realistic problem. 38 5.2.1 Chemotaxis phenomenon The study of chemotactic behavior of bacteria originated in the nineteenth century. In 1884, Pfeffer discovered the chemotaxis phenomenon in bacteria by observing a capillary tube which contained an attractant and motile bacteria. The bacteria accumulated near the higher concentration of the attractant. In 1901, Rothert, Jennings and Crosby found that bacteria often avoided the region of lower density in a capillary tube. In 1972, H. C. Berg and D. A. Brown analyzed the microscope motion of the bacteria by comparing the quantitative run-twiddle analysis between a wild type, a nonchemotactic mutant and an uncoordinated mutant. They discoverd that bacteria responsed to the external environment by modulating their pattern of motion shifting bewteen swimming and tumbling: swimming is a nearly straight run, tumbling is a very rapid and jerky movement which redirects the cell to a new random direction. Usually, the motion of swimming is followed by a brief tumbling. When bacteria encounter a favorable external environment, they decrease the tumbling probability and extend the time for swimmimg in the preferred directions. It results in a bias in the random walk and net movement in the gradient of the attractant. That is chemotaxis [8, 16]. The swimming and tumbling modes are modulated by the bacteria flagella. Each one consists a long filament, a left handed helix of protein subunits and a basal end attached to the cell membrane. The basal end serves as a turning wheel. It turns the filament as a propeller. The motion of swimming is produced by the anticlockwise rotation of the left-handed helix. It results in a forward movement. The clockwise rotation gives rise to a backward motion, producing tumbling. Bacteria sense the change of the external environment via membrane receptors 39 binding the chemical substrate. The cells reflect the change in the concentration of chemicals by the relative number of occupied binding sites, not by the absolute concentration of the substrate. “The receptor proteins perform as reporter molecules, which sending information about the current concentration of ligand to a central mechanism (tumble regular). The tumble regular then creates signal to rule the direction of the rotation of flagellar rotary motors” [8]. It is worthwhile to mention that we can derive the chemotaxis coefficient from the proportion between the occupied (denoted by No) and free receptors (denoted by N f). The ratio of the occupied receptors and the total number of receptors is NIIWXLW- Let T0 and Tf represent the mean time of a receptor occupation and the mean time when the receptor is free. Note that Tf is inversely proportional to the concentration of the substrate and T0 is a constant which is related to the process of internal cellular. Thus No To C ——=—=—— 5.5 Nf+N0 Tf+T0 K+C ( ) and a 1v0 _ K ac 55%, + N0) ‘ (K + C)2 823’ (5'6) where K :- (£72,2—) is a constant, and {If—f5: is known as the “receptor law”. The receptor law says that the chemotactic reponse vanishes when all of the receptors have been occupied [16]. We tested our numerical scheme by using Julius Adler’s series of work [4, 5, 6, 22] in the quantitative investigations into the chemotaxis. His work has inspired and 40 influenced researchers in the area of molecular biology and mathematical biology, especially in the field of pattern formation. We will briefly describe his experimental work and the results of one of his papers in 1966 before demonstrating the numerical results. 5.2.2 Experiment results Figure 5.2: Photographs showing bands of E. coli in a capillary tube. In all the ex- periments reported here, capillary tubes were filled with a liquid medium, inoculated at one end with 200000 to 2000000 bacteria, and then closed at the ends with plug of agar and clay, all according to procedure described in full elsewhere. The tube were incubated horizontally at 37C. The origin, which is turbid because of the bacteria that have not moved out, is visible at the left, then the second band of bacteria, then the first band. Plugs at the ends are not shown. The concentration of galactose was 0.00025 mole per liter. Julius Adler and his associates use E. Coli in most of his experiments. Since “ the vast knowledge of its biochemistry and genetics could be brought to bear on the problem.” It is an inhabitant in the human or animal intestine and it is nonpathogenic in normal environments. Ecoli is a. mobile bacteria, it moves by propelling itself with flagella around cells. E. Coli is known to be moving preferentially toward oxygen. 41 E. Coli were placed at one end of a closed capillary tube. The other end of the capillary tube was closed with plugs of agar and clay. It contained different amounts of oxygen and an energy source. Two clear and visible bands of bacteria with a constant speed moved from the end after adding bacteria into the tube. Some cells still remained at the end of the tube. If the amount of dissolved oxygen was lower than the energy source (the amount of oxygen is not enough to oxidize all the energy source), then the first band traveled along, consuming all of the oxygen. The second band would move out when the cells were able to deplete the energy source anerobically. On the other hand, when there is enough oxygen for the cells to oxidize the energy source, the first band totally exhausted the energy source aerobically. Under these circumstances, the second band used up the remaining oxygen to oxidize the endogenous energy source. In this situation, the bacteria in the first band generated a steeper gradient along the concentration of the energy source than the gradient in the concentration of oxygen. Bacteria then moved preferentially toward the higher concentration of the energy source. Bacteria in the second band created gradients in the concentration of oxygen and galactose. Bacteria therefore migrated toward the higher concentration of oxygen [4, 18, 24]. These two crowded regions of bacteria are the effect of chemotaxis. Since if random walk or diffusion is the only transport mechanism for the movement of bacteria, the distribution of bacteria density should look like a linear function which falls from high-to low density. 42 '8' 3‘ 8 l 1 Oxygen remaining (%) s ‘ 4 , a SECONd First Band Band Centimeters Figure 5.3: Utilization of oxygen by bands of bacteria in 0.00025 M galactose. Oxy- gen was measured polargraphically by inserting an oxygen needle electrode into the capliiary tube in 4-mrn steps when the first and second bands were visible where shown by arrows. 43 IOO - :5: .E’ 75 ~ .5 (U E '5 50 - a: O E a 25 '- (D l m l L L l c O I 4 I 8 Second First Band Band Centimeters Figure 5.4: Utilization of 0.00025 M galactose. The tube contained C14- galactose (0.00025 mole/ liter and 1000000 counts per minute per milliliter). At 4 or 8 hours, when the first and second bands were visible where shown by arrows, the tube was factionated into ten compartments, each 8 mm long. the contents of each compart- ment were chromatographed on paper, with n-butanol, acetic acid, and water (12:3:5) as the solvent, and the radioactivity in the galactose region of the chromatogram was measured in a paper-strip counter. 44 I00 I I l _ r T | l I | l l l l l 100 >- - :\°‘ ’3 75— - 0: 2i 5 75 i- . I» .E .S N E E cu so - E a, 50 l' .. 9 w c 2 a, U 2 a- - % x (D 25 ~ _ O I 0 f 5 4' '6 .T D 1 I l l l I 1 I 1 Second First 0 1 4 _' i a Band Band 590071" 3309 First Band Centimeters Centimeters Figure 5.5: Utilization of oxygen by bands of bacteria in 0.000005 M galactose. Mea- surements were made as in Figure 5.3 and Figure 5.4. Utilization of 0.000005 M galac- tose. The tube contained Cl4-galactose (0.000005 molde/ liter and 200000 counts per minute per milliliter). At 4.75 hours when the first and second bands were visible where shown by the arrows. 45 5.2.3 Mathematical Model for pattern formation in E.coli experiment We deal with a continuous, deterministic model to describe bacteria growth. The characteristics of the model we describe are all two dimensional. First, we will explain the following notation and its units H o u : bacteria density : [u]=cells/cm3. o g : concentration of galactose : [g]=mmole/cm3. o s : concentration of oxygen : [s]=mmole/cm3. o t : time : [t]=hours. o a: : distance along the tube (peri dish) : [x]=cm. o X1(g) = chemotactic sensitivity for galactose : [X1(s)]=cm2/hour. o X2(s) = chemotactic sensitivity for oxygen : [X2(g)]=cm2 / hour. 0 Du > 0 : bacteria motility : [Du]=cm2/hr 0 D9 > 0, and D3 > 0 : constant diffusivities : [D]: cm2 / hr. 0 Yg > 0 : bacteria yield constant, which gives the yield of bacteria per galactose taken up : [Y9]: mass u/mass g. 0 Y3 > 0 : bacteria yield constant, which gives the yield of bacteria per oxygen taken up : [Y3]: mass u/mass s. For a closed region, the time rate of change 'of the bacteria concentration should be equal to the bacteria flux across the boundary of the closed region. Assume that 46 reproduction is not negligible. This gives an 8"] CC ' 3) B—t— —— —% + BZT‘th , (5.7) where J = J (x, t) is the flux of the bacteria. It consists of two parts : diffusion and chemotaxis. The diflusion or random walk is given by —DuVu. (5.8) The chemotaxis flux is uX1(g)Vg + ux2(s)Vs, (5.9) where X1(g) and x2(s) are the chemotactic sensitivity to the galactose and oxygen respectively. They represent the strength of chemoattractant. The presence of u in the above equation is because the bacteria flux is the number of bacteria crossing a unit area per unit time. Therefore, increasing the number of bacteria will in- crease the bacteria flux. In the case of repellent, the flux will have a positive sign [16]. Here we use the Keller-Segel model with p=2 or the receptor law [18]. Co XKS(5) = m (5.10) where CO and C1 are positive constants. Note that as lz'ms —+ ooXKs(S) = 0, which means that bacteria will prefer to stay where they are if there is plenty of chemoattractant. 47 Thus, the total bacteria flux J is: J = uX1(g)Vg + ux2(s)Vs — DuVu. (5.11) For the growth term, we adapt Monod’s law V19 fl(g) = m: (5.12) where V1 > 0 and K1 > 0 are constants. If the amount of nutrient (galactose, in this case) is restricted, bacteria will consume as much as food as they can. When the supply of nutrient is abundant, the consumption rate is independent of the concentration of nutrient. We say the growth rate is saturated with respect to the nutrient [18]. Under this circumstance, the reproduction of bacteria only depends on the bacteria density and the eating rate per bacteria is Together with (5.11),(5.12), the bacteria equation u is ut = DuAu — V(uxl(g)Vg) — V(ux2(s)Vs) + f1(g)u, (5.13) where A = 82/63:? + 02/8xg is the standard Laplace operator, which simulates the diffusion process. V = 6/ 01:1 + 0/(9x2 is the gradient operator. The model also contains nutrient (galactose) and oxygen. Nutrient was taken up by bacteria at the same rate as bacteria reproduction. This suggests that the nutri- ent consumption (following Monod’s law) has a sigmoid characteristic. Meanwhile, 48 the nutrient also diffuses itself. A simular assumption also holds for oxygen. How- ever, oxygen does not contribute to the growth of bacteria. Hence, the equations for nutrient and chemotaxis are: gt = DgAg - f1(g)u/Yg, (5.14) st = DsAs -— f2(s)u/Y;. (5.15) where Yg > 0, Y; > 0 are the bacteria yield constants, Yg gives the yield of bacteria per nutrient taken up. Y3 gives the yield of bacteria per oxygen taken up. According to Adler’s experiment, the ends of the capillary tube are closed. No sub- strate flows through the boundary. Therefore, the normal derivative on the boundary is zero. i.e, 811 89 83 anlan an an anlan 0 (5 6) Moreover, the initial condition for bacteria is : uo(:z:) if “an” S K u(:1:, 0) = 0 Otherwise. which approximates the initial inoculum of the bacteria in the center of the agar plate (in the two dimensional case). Galactose and oxygen are uniformly distributed at the beginning. The initial condi- tions fitting Adler’s experiment are 9(33, 0) = 90(1), (5.17) 49 5(23, 0) = 30(2). (5.18) 5.2.4 The coefficient function To complete the model equations, we have to estimate the parameter function. It is crucial to have an accurate parameter estimation to ensure the presentation of patterns. The original paper does not contain the quantitative data of the experiment, so we use ChiChia Chiu’s paper [14] as a reference to guess the correct order for the parameter function. The values of parameters are chosen as follows. 0 Du = 0.001, D9 = 0.005 D3 = 0.033 0 u0(x) = 5 x 10‘6 if ”2:” S 6, where 6’ = 7.5 x 10‘2 0 30(2) 2 4 x 10‘5, 90(2) 2 2.64 x 10‘5 . Y9 = 0.001, Y3 = 0.0001 CKI = 0.02,DK1= 2 x 10—6 CK2 = 0.04, DK2 = 3.3 X 10—5 v1 = 0.35, K1 = 4 x 10-6 V2 = 0.6, K1 = 5.5 x 10-5 Q = [—R, R] x [—R, R], where R =1 These parameters come from experiments and literature. In the numerical sim- ulation, we fix all parameters except u0(:r), 90(3) in order to accommodate Adler’s experiment. 50 5.2.5 Intuitive explanation of pattern formation mechanism The nutrient (galactose) and chemoattractant (galactose, oxygen) will exhaust as time goes to infinity. This implies that the growth of bacteria and the effiect of chemotaxis will eventually become zeros. In the end, the only term which remains in effect is the diffusion of bacteria. The solution become spatially homogeneous. This is not an interesting phenomenon in pattern formation. Thus, we have to investigate a dynamic solution [20]. Note that galactose not only plays a role as an energy source but also as chemoattractant. Therefore, galactose and oxygen both serve as chemoattractant. They both direct the movement of bacteria. However, oxygen does not contribute to the growth of bacteria. The bacteria take up suger and oxygen, creating a gradient of nutrient and chemoattractant. They move preferentially toward higher concentration of substrates and generate a population there. However, the diffusion of the bacteria has a dissipated effect, which flattens out the population peak. The speed of diffusion and the eating speed of the bactera for these chemicals are not identical. As a consequence, two sharp and visible bands migrate out. 5.2.6 An optimal two-stage numerical scheme Let Q = [0,8]2 with mesh size h = 26/ N. Let $ij E Q be a grid point, then $i+1,1 = $21 + h and $1,j+1 = $1j + h, etc. The time step is denoted as At, and “ij % “$.5- Let us solve the equations for the nutrient concentration, 9, the attractant con- 51 Ill-i. “-1. "I! . centration, s, using the optimal two stage scheme: 89 92'+1,' —' 293' + 91—1; 92', '+1 — 2923' + 92', '—1 =D,( , '3 .)+,,g(. , ,) 797 h? h? - f2(gijuij)/Y:q D9 D9 4Dg D9 D9 = figz‘uj + 72—29041 — 7923' + 32-92;” + 3590—1 - f2(gz'juz'j)/Yg' (519) 33 5i+1,j - 28m + 3i—1,j 52',j+1 — 28M + 82332—1 a : D3( h2 ) + D3( hz ) - f1(50'uz'j)/Ys D3 D3 4D3 D3 D8 = 7133i+1j + 72792741 — 7513' + fist—U + fisz’j—l " f1(3ijuij)/Ys- (5.20) The boundary conditions are 9—1j = 913', 9N+1j = gN—lja gi—l = 921, and 9—1,N+1 = 9i,N—1, 8—1j = 813', 3N+1j = «SN—13', Si—l = Sn, and 8—1,N+1 = 5i,N—1a where i,j = 0,1,2... ,N. Now, we can investigate the time step of the ODE system 8./8t = M (t) + f (:r, t) for g and s. From (5.19) and the boundary, the diagonal entries of 9 matrix and the summand of the off diagonal entries are h2 ’ h2 52 respectively. By the Gerschgorin Theorem, we get 40, 4D,, _ 757' [A+ )1? Is For 2' = 0, 1, 2, - - - ,N. Therefore, the real part of the extreme eigenvalues of M for g are contained in —8D [ hg ,0]. Hence, the optimal time step for g is h2 = —, .21 4,,9 (5 > Atg Similarly, the optimal time step for s-equation is h2 At, = 40 . (5.22) In the following section, we present optimal two—stage scheme for solving the reaction-diffusion-chemotaxis equation. The equation for it has the following form at = DuAu — V(uxl(s)Vs) — V(ux2(g)Vg) + f1(g)u. (5.23) The optimal two—stage scheme for u is defined by the following: Bu _ Amth-j Ayhuij at _ D" h? +D“ h? X1(3)V$s+lxl(s)sz| _ ( 2h Vzrh)ijuij ”(3)sz — [X1(5)V$s| - ( 2h Vihlz‘juij (5.24) 53 X1(8)\7 8+ IX (8W 8| _ ( y 2h 1 y Vyhlz'juz'j X (8W 8 - IX (8W 8| ( 1 y 2h 1 y Vyh)z.7u 3.7 X 9V59+ X ngg _ ( 2( ) 2h] 2( ) IV mh)ijuij X2(9)V x9_ IX2(9) ngl ( 2h th)iJu ui] x (9)Vg+lx (9)V 9| .. ( 2 y 2h 2 y Vyh)ijuij X2(9)Vy9_lX2(9 )V ygl+ ( 2h V3 0sz 0 + fl(gij)uij (5.25) where vgh and V3}; are first order backward difference operators such that (Vii. + Vgh)u5= (Uij_ui—1j)+(uij - uz‘j—l), and V: and V+ h 3/ h are first order forward difference operators such that (Vi); + vyh)uij= (ui+1j— Uij) + (“'in — Uij) V2,, Vy are the centered difference operators Va: = ui-i-lj — ”Mi—13', Vy = uij+1 — uz'j—L Note ._ + - Axhuz'j — vxhuij — Vrhuij' and (X1(3))ij = X1(3ij)- The zero flux boundary conditions for u are approximated by 54 71—13“ = U133 UN+1j = uN—lja ui—l = mi, and u—1,N+1 = ”LN—1' Therefore, (5.24) can be written as % ___ DuUi+1j — 2:? + Uni—13’ + Baum-+1 — 2:3]: + 11.25-} _ (X1(S)V$s :hIX1(s)V$s| )ijwij _ iii—13') _ (xl(s)v.s 215' X1(s)v,s| )ij (qu _ W) — (X1(S)Vys :Ile1(S)VySI lifluzj — ”(j—1) — (X1(8)vys 2th1VySI )ij(uz'j+1 — uz’j) _ (X2(9)ng :th2(9)ngl)ij(uij _ qu) _ (X2(9)ng ;h|X2(9)ng| )ij(ui+1j _ 0.3-) _ (X2(9)Vy9 :hIX2(9)Vy9l)ij(uij _ Haj—1) _ (X2(9)Vyg ghl><2(g)Vy9|>150,”+1 _ W) + f1(gij)uz'j- That is, _8_u = [9—1, + (X1(s)sz + IX1(3)Vx3|)H at h2 211 ‘3 x ng+ x ngg ( 2() 2]] 2() 55125—1.- fi (X1(5)Vy5 + |X1(S)Vy3l)__ h? 2h ’3 X2(9)Vy9 + IX2(9)Vy9| .. ( 2h )0] ufi-l —4Du 20(3)sz + lX1(S)vx3l __ —h—2_ — ( 2h )” + l 55 X1($)Vy3 + |X1(3)Vy3| - ( + (X1(3)Va:8 — |X1(5)Va:5l 2h 2h Du 2h ) 2h ) _ (X2(9)Va:9 + [X2(9)Va:9l) (X2(9)Vzg - IX2(9)V59|) ”(3)sz — |X1(s)sz| (X1(8)Vy8 — |X1(3)Vy5|),_ 2h ’9 _ (X2(9)Vy9 + |X2(9)Vy9|) 2h X2(9)Vy9 "' lX2(9)Vygl + ( 2h )0] 1.7 U + [75 _ ( 2h (>901)ng - [X2(9)Va:9| 2h D )z'j )z'jluz'm X1(3)Vy5 — lX1(3)Vy3l + [77.7 _( 2h (X2(9)Vy9 — IX2(9)Vygl 2h where i,j = 0,1,2~- ,N. )z'j )ijluz'j+la Let (0)” denote the time level. Then the diagonal entries of the spatial matrix M .= h2 _( = M(t") +2M(t"+1) —4Du X1(8)Va=8 + [X1(3)Vx3| of the ODE system for u is (X1(5)V:c5 + |X1(3)Va:3l)n+1 4h X1(5)Vx5 - |X1(5)Vx5| 4h X1(s)Vys + IX1(s)Vys| 4h X1(5)Vy3 - le(3)Vy5l 4h _ (”WWW + |X2(9)ngl 1... 4h '3 (X2(9)ng - [X2(9)ngl )1. 4h 2’ _ (X2(g)vyg + lX2(g)V3/gl)n 4h (X2(g)Vyg - |X2(9)Vy9| )9 4h 2] +( -( + ( )2;- + )1} — )2; + _ (20(9)va + lX2(g)V:cgl)n+l +( fj— +( )3 _ 4h ‘7' 20(3)sz — lX1(S)vI€Sl n+1 X1(3)Vy3 + lX1(3)Vy3l n+1 (newts — lxl(s)Vysl ...1 4h z] 4h 27' 20(9)ng - lX2(9)Vx9l)n+1 4h '7' (X2(g)Vy9 + IX2(9)Vy9| )1.“ 4h '3 20(9)va _ lX2(g)Vy9l)n+1 4h ’9 ° The summand of the absolute value of the off diagonal entries for each row is Z laijl j_laj#i 4D,, h2 -( 20(3)sz "' lX1(S)Va:5l 4h 20(5)sz — |X1(5)vxsl - ( 4h X1(3)Vx8 + |X1(S)Vz3| + ( 4h X1(3)Vy5 - |X1(8)Vy8| - ( 4h X1(3)Vy3 + |X1(3)Vy3l + ( 4h )gfl >:;-+ )3- ) —( +( -( + X2(g)v:cg - lX2(g)V:cgl)n, __ 2] X2(g)vrcg + lX2(g)vxgl)n +( 4h 4h X2(9)Vyg — IX2(9)Vy9| (X2(9) 4h V319 + [X2(9)Vy9l 4h n U 2'1 )3. )3} -( +( >009)sz + le(s)sz| ( 4h X1(3)Vy3 _ lX1(S)V3/Sl )n+l )0.“ 1] 4h X1(5)Vy3 + lX1(3)Vy3l ii ( -( +( 4h )flfl 2] X2(9)V29 - IX2(9)V:::9| )7.“ 4h X2(9)V:cg + lX2(g)V$gl )n+l ij 4h 27’ X2(g)vyg — lX2(g)V3/gl )n+l 4h 20(9)va + lX2(g)V3/gl )n+l ij 4h ij' Therefore, the estimated largest and the smallest eigenvalues can be found by the following. The smallest eigenvalue of the spatial matrix M is “_amin : ( ( ( ( ( —8Du h2 X1(5)V23 + |X1(3)V53| -( X1(s)sz + IX1(s)V$s| 2h X1(3)Vy3 + lX1(3)Vy5l 2h X1(3)sz - |X1(3)V:c3| 2h ) X1(8)Vy3 — |X1(3)Vy5| 2h 2h X2(9)ng + |X2(9)V229l 2h 57 )n+1 2'J' n — ij >23+ >:3-+ )2;— ) Tl 1'} X1(8)Vy3 + |X1(3)Vy5| ( 2h X1(5)V15 - le(3)Vx3l ) ( 2h X1(3)Vy3 - |X1(S)Vy8| ( ( 2h ) X2(9)ng + |X2(9)ngl n + 1 “ij )0.“ U n+1 ij )n+1 2h ij _ (X2(9)Vy9+|X2(9)Vy9l)n._ (X2(9)Vy9+lX2(9)V1/9l 541 2h ‘9 2h ‘3 X Vx—X V597. gixg—gixgn + ( 2(9) 92hl 2(9) |)ij+( 2( ) 2h] 2( ) It)“ X2(9)Vy9 - IX2(9)Vy9| n X2(9)Vy9— [X2(9)Vy9l n+1 + Note that “amin < 0. The largest eigenvalue of the spatial matrix M is —amax = 0. Therefore, the time step for b system is 4D,, lsts+ lsts n M) = 1/(7+(X( ) 4hlx () I)” 20(3)sz + lX1(5)V$3l)n+l 4h ‘7 X1(3)Vy3 + lX1(3)Vy5l )1}. (X1(5)Vy3 + lX1(3)V3/Sl )n.+1 4h ’3 4h '3 _ (>009)sz — |X1(3)Vx3[ )5 _ (>00)sz - |X1(3)V:c3| )1...“ 4h ’3 4h ‘3 _ (X1(8)Vy8 — |X1(8)Vysl)n. _ (X1(3)Vy3 _ lX1(8)Vy'sl)T-r+l 4h ‘3 4h ’3 (X2(9)ng + |X2(9)ngl )7. + (X2(9)V59 + |X2(9)ng| )7.“ 4h. ‘9 4h '3 _,_ (X2(9)Vy9 + IX2(9)Vy9|).._ + (X2(9)Vy9 + [X2(9)Vy9|)p_+1 4h “J 4h '3 (X2(9)ng - IX2(9)V59I),. _ (29(9)ng - IX2(9)V59|)..11 4h '3 4h '3 _ (X2(9)Vy9 - IX2(9)Vy9| )5 _ (X2(9)Vy9 - |X2(9)Vy9l)15+1 4h 23 4h '1 ' + ( + ( The optimal time step for the our model equation is min {Atb, Atg, Ats} . (5.26) 58 5.2.7 Results of numerical experiment In this section, the results of the numerical experiments will be demonstrated to compare with Adler’s experimental results [4] This picture is taken from Adler’s experiment. The explanation of this picture will First Band ‘ . ll 0 fir“ -< Centimeters at 4 hours Second Band N0 second band .L_J o 0.5 LO fi' 210 5.0 0.001 M Galactose Figure 5.6: The relation between the amount of galactose and the distance of the bands of bacteria. be described upon each graph of numerical simulation. We will use 2D pictures to simulate the 1D phenomenon. The reason is, first, Adler’s experiment can be done in the agar plate. Secondly, most of patterns are shown by using 2D pictures. We want to make sure that our numerical method is suitable for general patterns. 59 -7 x 10 1 8 0.5 6 0 ‘ ‘ 4 I —0.5 2 -1 0 —1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Concentration of Galatcose is 0.00025 Concentration of galactose is 000046 O -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Concertration 01 Galactose Is 0.0008 Concentration of Galactose is 0.001 Figure 5.7: This set of pictures simulates the movement of bands of bacteria when the concentration of galactose is higher than oxygen. The amount of galactose and oxygen varied during the experiment. When the concentration of galactose was above a certain level, the second band moved less far from the origin as the concentration of galactose was increased. It showed that the band would not proceed unless the bacteria ingested all the galactose. On the other hand, the distance for the first band did not vary. Probably, the consumption rate of oxygen was saturated. 60 -7 x210 1 1.8 0.5 1.6 1.4 0 1.2 -0.5 1 _1 0.8 -1 0 1 -1 -0.5 0 0.5 1 Concentration of galactose is 0.00004 Concentration of galactose is 0.00005 -7 x310 “ 2.5 2 1.5 1 -1 0 1 -1 0 Concentration of galactose is 000006 Concentration of galactose is 0.000099 Figure 5.8: The picture shows the movement of bacteria when the concentration of galactose is below the oxygen. The experiment also tested the movement of bacteria when the concentration of galactose was under a certain amount. In this case, increasing the concentration of galactose slowed down the movement of the first band and facilitated the movement of the second band. Since bacteria in the first traveled along and consumed all galactose, it would take more time for the bacteria to consume nutrient when the amount of galactose was increased. At the same time, the consumption of oxgen was also increased for the first band of bacteria. Therefore, less oxygen was left for the bacteria in second band. Bacteria in the second band would then needed less time to use it up. 61 O ’-1 05 o 05 1 '-1 . . -0.5 0 0.5 Concentration 01 oxygen is 0.001 Concentration of oxygen is 0.008 Figure 5.9: Relation of the amount of oxygen and the distance of the bands of bacteria. When the concentration of galactose was too high, only the first band appeared during the experiment. The experiment of different oxygen levels was tested under this condition. It was observed that the higher the concentration of oxygen, the slower the band moved. The condition was similar to the galactose, the bacteria moved out only when they used up all the surrounding oxygen. Therefore, increasing the concentration of oxygen would decrease the speed of the band. 62 CHAPTER 6 Conclusion This thesis develops a new numerical method to solve the reaction-diffusion- chemotaxis systems. This optimal two-stage method is to improve the performance of the explicit two-stage Runge-Kutta method for the time evolution problems. It retains the advantage of the explicit two-stage Runge-Kutta method of simple, easy implementation, and higher accuracy. In addition, the explicit formulation of the time step can be calculated. The time step is adapted and locally optimal. The consistency, convergence and the stability issues have been proved. The comparison between the optimal two—stage scheme and the two-stage Runge-Kutta method has been made. Table 1 and Figure 1 show that the optimal two-stage scheme is accurate and converges faster than the RK method for the time evolution problems. A real problem is also used to test out numerical method. Mathematical equations are con- structed to model Adler’s experiment in which patterns of two bands emerge. The numerical simulation and the experiment results are well matched. 63 BIBLIOGRAPHY 64 BIBLIOGRAPHY [1] D. A. Anderson, J. C. Tannehill, R. H. Pletcher, Computational fluid mechanics and heat transfer, Series in computational fluid mechanics and heat transfer, Hemisphere Publishing Corporation. [2] VanBuskirk AM, Burlingham WJ, J ankowska-Can E, Chin T, Kusaka S, Geissler F, Pelletier RP, Orosz, Human allograft acceptance is associated with immune regulation, J Clin Invest, Vol 106, 145-155, 2000. [3] U. Ascher, S. Ruuth, R. Spiteri, Implicit-explicit Runge-Kutta methods for time- dependent partial differential equations. Special issue on time integration. Appl. Numer. Math., 25, 151-167, 1997. [4] Julius Adler, Chemotaxis in Bacteria, Science, New Series, Vol 153, No.3737. Aug. 12, 1966, pp.708-716. [5] J. Adler and B. Templeton, The eflect of environmental conditions on the motility of Escherichia coli, J. Gen. Microbiol. 46:175—184, 1967. [6] O. Baracchini and J. Sherris, J. Pathol.Bacteriol, 77, 565, 1959. [7] E. O. Budrene, H. C. Berg, Complex patterns formed by motile cells of Es- cherichia coli, Nature, 349, 630 - 633, 1991. [8] Biology of the C'hemotactic Response,edited by J .M. Lackiie, Cambridge Univer- sity Press, Cambridge, 1986. [9] H. C. Berg, Random Walks in Biology, Princeton University Press, Princeton, NJ, 1989. [10] E. O. Budriené, A. A. Polezhaev, and M. O. Ptitsyn, Mathematical modeling of intercellular regulation causing the formation of spatial structures in bacterial colonies, J. Theor. Biol. 135, 323-341, 1988. 65 [11] Chichia Chiu, Optimal one-stage and two-stage schemes for steady state solutions of hyperbolic equations, Applied Numerical Mathematics, 11, 1993, 475-496. [12] Mei-Qin Chen, Chichia Chiu, Optimal m-stage Runge-Kutta Schemes for steady- State Solutions of hyperbolic equations, Numerical methods for partial differential equations, 9, 643-666, 1993. [13] Chichia Chiu, David A. Kopriva, An optimal Runge-Kutta method for steady state solutions of hyperbolic systems, SIAM J .NUMER. ANAl, Vol 29, N02, pp 425-438, 1992. [14] Chichia Chiu, Frank Hoppensteadt, Mathematical models and simulations of bacterial growth and chemotaxis in diflusion gradient chamber. J. Math. Biol, Spring-Verlag, 2000. [15] G. W. Gear, Numerical initial value problems in Ordinary Difierential Equations, Prentice-Hall, Inc. 1971. [16] Eshel Ben-Jacob, Inon Cohen, Ido Golding, Yonathan Kozlovsky,Modeling branching and chiral colonial patterning of lubricating bacteria, Mathematical models for biological pattern formation : Frontiers in Application of mathemat- ics, edited by Philip K. Maini, Hans G. Othmer, The IMA Volumes in Mathe- matics and its applications Vol. 121, Spring-Verlag, 2001. [17] P. Grindrod, The theory and applications of reaction-diffusion equations: pat- terns and waves, 2nd Ed, Oxford University Press Inc., New York, 1996. [18] E. Keller, L. Segel, Model for chemotaxis, Jour. Theor. Biol.,30, 225-234, 1971a. [19] K. W. Morton, Numerical Solution of convection-diflusion problems, Applied Mathematics and Mathematical Computation 12, Chapman & Hall. [20] J. D. Murray, Mathematical Biology. 3rd edition, Springer—Verlag. [21] M. R. Myerscough, P. K. Maini, K. J. Painter, Pattern formation in a generalized chemotactic model, Bulletin of Mathematical Biology, 60, 1-26, 1998. [22] J. C. Sherris, N. W. Preston, J. G. Shoesmith, J. Gen. Microbiol. 16, 86, 1957. [23] G. D. Smith, Numerical solution of partial differential equations: finite diflerence methods - 3rd Ed. Oxford applied mathematics and computing solutions. [24] Scribner, T., Segel, L., Rogers, E.,A numerical study of the formation and prop- agation of traveling bands of chemotactic bacteria, J. Theor. Biol, 189-219, 1974. 66 [25] Thring, A., The chemical basis of morphogenesis. Proc. R. Soc. Lond. 237: 37-72, 1952. [26] Rebecca Tyson, L. G. Stern, Randall J. LeVeque, Fractional step methods applied to a chemotaxis model, J. Math. Biol. 41,455—475, 2000. [27] R. C. Tyson, S. Lubkin, J. Murray, Model and analysis of chemotactic bacterial patterns in a liquid medium, J. Math. Biol, 38, 359-375, 1999. 67