--\.’ I‘v _ IIIllfillljlflllillllfillilllllll 1686 2371 This is to certify that the dissertation entitled Engineering Applications of Microbial Chemotaxis presented by Mark Thomas Widman has been accepted towards fulfillment of the requirements for Ph.D. degree in Chemical Engineering K 7715M Zb/fl/Oék Major professor S/HK/C/7 Date MSU is an A/lirnmtiw- At'ttnra/It'qiml Opportunity Institution 0-12771 LIBHMY Mlchlgan State Unlversity PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE 1/98 cJCIRC/DateDue.p65—p.14 ENGINEERING APPLICATIONS OF MICROBIAL CHEMOTAXIS By Mark Thomas Widman A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Chemical Engineering 1997 ENG Chemo‘ chemicals call chemotaxis to to develop or r tnolsand meth the laser dit’t‘ra Sacral mathet experimental o predictions we modeling simt Chttmotacttc rt saturation calcl The Set Which Chemotz fully take adv; all)Plicattong th 3”" bi(tremedi ABSTRACT ENGINEERING APPLICATIONS OF MICROBIAL CHEMOTAXIS By Mark Thomas Widman Chemotaxis is the ability of an organism to move in response to gradients of chemicals called chemoattractants. Many species of bacteria are known to exhibit chemotaxis to a wide variety of chemoattractants. The first objective of this research was to develop or refine tools and methods to study the chemotactic response of bacteria. The tools and methods include the diffusion gradient chamber (DGC) and its associated units, the laser diffraction capillary assay (LDCA), and microelectrodes and microbiosensors. Several mathematical models of the DGC system were written that allowed predictions of experimental outcomes to be made without running the actual experiments. The modeling predictions were verified by comparison to experimental data. Methods to analyze the modeling simulations were developed to further enhance the understanding of the chemotactic response. These analysis methods included bacterial flux and receptor saturation calculations. The second objective of the research was to identify applications or systems in which chemotaxis was beneficial, and then to develop ways to engineer those systems to fully take advantage of the benefits offered through the chemotactic response. The three applications that were selected were microbial competition, selection of mutants, and in situ bioremediation. Chemotaxis was predicted by the mathematical model to impart a l «why—".1 , '1’” _-._ -fi'a- competitive ac' preliminan' W testunse wt (alt. In situ themotactit t {initiation of telngused tor chemotaxis to response of P KC and ESt‘iit’ competitive advantage to a bacterial population in certain non-mixed environments, and preliminary work was done to experimentally validate the predictions. The chemotactic response was used as a selection agent in the isolation of desirable mutants of Escherichia coli. In situ mutagenesis was performed on feedback inhibited E. coli, and the chemotactic response was used to separate feedback resistant mutants from the population of inhibited bacteria. Pseudomonas stutzeri strain KC, a bacterium currently being used to remediate a carbon tetrachloride contaminated aquifer, was found to exhibit chemotaxis to chemicals present in the aquifer. Experimental protocols to study the response of Pseudomonas KC around objects, and competition between Pseudomonas KC and Escherichia coli were developed. The at tori reported performed in; tort on the responsible to 1 \soul with me on 2 those researc Marshall Bre( Dr. .\1 {hate reache DGC and me microorganisr ADI algorithr and Faith Pep IOWe thank my wif has Worked e: her devotion r ACKNOWLEDGMENTS The efforts of many people have been instrumental in helping me to complete the work reported in this thesis. I am grateful for the help of a number of undergraduates who performed many of the experiments. Laura Booms and Petty Setiawan both did good work on the Pseudomonas KC swarm plate experiments. Sebastian Schmidt was responsible for making the Laser Diffraction Capillary Assay a reality. I would like to thank my fellow graduate students, Tyler Ames, who collaborated with me on a project that impacted both of our research programs, and Mark Mikola, whose research on mutant selection led to another use for chemotaxis. Both of them and Marshall Bredwell provided me with many insights into my research. Dr. Mark Worden has mentored me through my four years, and with his guidance I have reached the goals that I set for myself. Dave Emerson taught me how to use the DGC, and more importantly, gave me a good introduction to handling and thinking about microorganisms, from a non-engineering viewpoint. Chichia Chiu supplied me with the ADI algorithm used to solve some of the math models. Candy McMaster, Julie Caywood and Faith Peterson in the graduate office have also been invaluable resources and friends. I owe a lot to my Mom and Dad, for everything. Most importantly, I would like to thank my wife, Shelly. She has put up with me all through my graduate school years, and has worked extremely hard so that I could stay in school. I hope to finally be able to repay her devotion to me by proving to her that all this was really worth it. iv LIST OF TAIE LIST OF FIG LIST OF NO.‘ 1. L\TRODL" 3. OBJECT 1\ 3.3 DevelO' 3.3 Appliez 3. CHESIOT: 3.1 The tilt 3.3 Experir 4. DEVELOP 4.1 Improv 4.1.1 1m. 4.1.2 Qu 4.1.3 Me 4.1.4 Fag 4.1.5 Re 4.1.6 Mi 4.2 Derek) 4.2.1 B0 4.22 lni 4.2.3 Ex 4.2.4 C0 4.2.5 Mr 4.2.6 Re 4.2.7 Ef: 4.2.8 CC 4.2.9 A5 4-3Ana1yg 43.1 Fit 43-2 Re 4-3.4 CI. 4'4 Variati '5 Laserc 5' MICROBI TABLE OF CONTENTS LIST OF TABLES ............................................................................................................ vii LIST OF FIGURES .......................................................................................................... viii LIST OF NOMENCLATURE ........................................................................................... xi 1. INTRODUCTION ............................................................................................................ 1 2. OBJECTIVES AND SIGNIFICANCE ............................................................................ 3 2.1 Deve10pment of tools ................................................................................................. 3 2.2 Applications for chemotaxis ...................................................................................... 4 3. CI-IEMOTAXIS BACKGROUND .................................................................................. 6 3.1 The chemotactic mechanism ...................................................................................... 6 3.2 Experimental systems for measuring chemotaxis ...................................................... 7 4. DEVELOPMENT OF TOOLS ...................................................................................... 11 4.1 Improvements to the DGC method .......................................................................... 11 4.1.1 Image acquisition system .................................................................................. 11 4.1.2 Quantification of cells by grayscale analysis .................................................... 12 4.1.3 Method to measure membrane mass transfer coefficient .................................. 15 4.1.4 Faster approach to steady-state gradients .......................................................... 17 4.1.5 Recycle of source and sink flasks ..................................................................... 20 4.1.6 Microsensors and microbiosensors ................................................................... 20 4.2 Development of model ............................................................................................. 24 4.2.1 Boundary conditions ......................................................................................... 28 4.2.2 Initial conditions ............................................................................................... 30 4.2.3 Experiments used to validate model ................................................................. 30 4.2.4 Computer simulations ....................................................................................... 33 4.2.5 Modeling parameters ......................................................................................... 33 4.2.6 Results of validation simulations ...................................................................... 36 4.2.7 Effects of grid-spacing on simulations .............................................................. 40 4.2.8 Concentration profiles of chemoattractants and nutrient .................................. 40 4.2.9 Assumptions and simplifications in the mathematical model ........................... 47 4.3 Analysis of modeling results .................................................................................... 50 4.3.1 Flux calculations ............................................................................................... 50 4.3.2 Receptor saturation ........................................................................................... 55 4.3.3 Pattern dynamics ............................................................................................... 57 4.3.4 Chemotactic wave in response to oxygen ......................................................... 58 4.4 Variations on original model ................................................................................... 5 8 4.5 Laser diffraction capillary assay .............................................................................. 61 5. MICROBIAL COMPETITION ..................................................................................... 66 5.1 Compe‘ 5.3 Compe' 5.3 Modeli. 5.4 Analys: 5.4.1 Th1 5.4.3 Sat 5 5 Compe 5.5.1 E11 5.5.3 \‘a 5.5.3 \‘a 5.5.4 To 5.5.5 \'a 5.5.6 \‘a 5.5.7 Va 56 Future 6. SELECTK 6.1 Expert 6.3Mode1: BIOREME 7.1 Experi 7.1.1 Sv 7.1.3 D( 7.1.3 M. 7-3 Future 8. SL'MMAF APPENDIX APPENDIX APPENDIX APPENDIX APPENDIX APPENDIX APPENDIX LIST OF RE 5.1 Competition background .......................................................................................... 66 5.2 Competition mathematical model ............................................................................ 69 5.3 Modeling parameters ............................................................................................... 71 5.4 Analysis of competition simulations ........................................................................ 73 5.4.1 The dynamic competition factor ....................................................................... 73 5.4.2 The average concentration curve ....................................................................... 75 5.4.3 Saturation of chemoreceptors ............................................................................ 77 5.5 Competition simulations .......................................................................................... 77 5.5.1 Effect of modeling parameters on competition results ..................................... 77 5.5.2 Variation of ans ................................................................................................ 78 5.5.3 Variation of vas .................................................................................................. 82 5.5.4 Variation of KDAS .............................................................................................. 84 5.5.5 Variation of chemoattractant concentration ...................................................... 87 5.5.6 Variation of inoculation conditions ................................................................... 90 5.5.7 Variation of nutrient initial conditions .............................................................. 92 5.6 Future work on competition ..................................................................................... 95 6. SELECTION OF MUTANTS ........................................................................................ 99 6.1 Experimental work ................................................................................................... 99 6.2 Modeling selection of mutants ............................................................................... 100 7. BIOREMEDIATION ................................................................................................... 1 13 7.1 Experimental results .............................................................................................. 115 7.1.1 Swarm plates ................................................................................................... 115 7.1.2 DGC experiments ............................................................................................ 119 7.1.3 Movement of wave around obstacles .............................................................. 123 7 .2 Future work on bioremediation .............................................................................. 126 8. SUMMARY ................................................................................................................. 129 APPENDIX A Miscellaneous chemotaxis experiments .................................................. 131 APPENDDI B The Alternating-Direction Implicit Method ............................................ 135 APPENDIX C Instructions for FORTRAN model .......................................................... 137 APPENDIX D The FORTRAN model ............................................................................ 141 APPENDIX E Input files for model. ................................................................................ 184 APPENDIX F Matlab program files ................................................................................ 187 APPENDIX G Dimensionless model .............................................................................. 196 LIST OF REFERENCES ................................................................................................. 199 vi Table 1.Dilut Table 3. Para Table 3. Matt Table 4. Base Table 5. Base Table 6. Para LIST OF TABLES Table 1. Dilution factor and OD for grayscale calibration ................................................. 13 Table 2. Parameter values for validation of model ............................................................ 36 Table 3. Maximum predicted cell fluxes for Run 1 at 20 hours. ....................................... 51 Table 4. Base case parameters for competition simulations. ............................................. 72 Table 5. Base case initial conditions for competition simulations ..................................... 87 Table 6. Parameter values for inhibition model ............................................................... 106 vii i _" “.— ’Té'TtfiT Figure 1. Rar. Figure 3. Rut Figure 3. Dia Figure 4. .\'or Figure 5. Cal Figure 6. Det Figure 7. Th1 Frgure 8. Sin Figure 9. Th1 figure 10. E, Figure 11. L‘ Fitrue 12. st Figure 13. 31 Figure 14. 81 Figure 15. c Figure 16. T Figure 17 It Figure 18.0 Figure 19. G Figure 20. C LIST OF FIGURES Figure 1. Random walk ........................................................................................................ 6 Figure 2. Runs and tumbles in the presence of a gradient. .................................................. 7 Figure 3. Diagram of the diffusion gradient chamber. ......................................................... 9 Figure 4. Normalized grayscale vs. optical density for E. coli. ......................................... 14 Figure 5. Calibration curve for E. coli in DGC .................................................................. 15 Figure 6. Determination of the membrane mass transfer coefficient. ................................ 16 Figure 7. Three-slab pouring method. ................................................................................ 19 Figure 8. Simulation of approach to steady-state. .............................................................. 19 Figure 9. Three dimensional glucose microsensor calibration plane. ................................ 22 Figure 10. E. coli forming ring in response to glucose. ..................................................... 23 Figure 11. Use of microsensors, microbiosensors, and image analysis, DGC modeling...24 Figure 12. Schematic diagram of DGC, showing modeling coordinates. .......................... 32 Figure 13. Simulation compared to experiment for 1.0 mM aspartate source ................... 38 Figure 14. Simulation compared to experiment for 0.1 mM aspartate source. .................. 39 Figure 15. Cell concentration profiles across the two central axes of the DGC. ............... 40 Figure 16. Top: Typical cell profile. Bottom: nutrient profile. .......................................... 42 Figure 17 Top: Typical S profile. Bottom: Q profile. ........................................................ 43 Figure 18. Oxygen gradients in the DGC arena ................................................................. 44 Figure 19. Glucose gradients at several positions in the DGC .......................................... 46 Figure 20. Comparison of model with and without hyperbolic tangent term. ................... 49 viii 3":_.*—"*‘,,. m‘ Figure 31. F1 Figure 33. F1 Figure 33. F1 Figure 34. It Figure 35 To Figure 36. 52 Figure 37. F1 Figure 38. 11 Figure 39. C- Figure 30. I Figure 31 V; Figure 31 0 Figure 33_ v Figure 34. V Figure 35: v Figure 36. p. Figure 37, p Figure 38. D Figure 39. C Figure 40. C Figure 41_ C Figure 21. Fluxes due to response to S chemoattractant. ................................................... 52 Figure 22. Fluxes due to response to Q chemoattractant. .................................................. 53 Figure 23. Fluxes due to random motility. ......................................................................... 54 Figure 24. Total mass fluxes. ............................................................................................. 54 Figure 25 Top: 438 for Run 1. Bottom: dig for Run 2. Note different vertical scale ............. 56 Figure 26. Sample simulations from second-generation model. ........................................ 61 Figure 27. Filling the capillary in the LDCA ..................................................................... 63 Figure 28. The random motility coefficient of Pseudomas KC. ........................................ 65 Figure 29. Coordinate axes used for competition analysis. ............................................... 74 Figure 30. Typical simulation, illustrating average concentration curve calculation. ....... 76 Figure 31 Variation of XOAS- ............................................................................................... 79 Figure 32. Other results of variation of 960/43 ...................................................................... 80 Figure 33. Variation of vAs. ............................................................................................... 83 Figure 34. Variation of KDAS. Time=30 h for all figures. .. ................................................ 85 Figure 35: Variation of S source concentration. Time=28 h for all figures. ...................... 89 Figure 36. P0p. A inoculated farther from S source than pop. B. ...................................... 90 Figure 37. P0p. A inoculated to the side and away from S source ..................................... 91 Figure 38. Dynamic competition factors for two inoculation patterns. ............................. 92 Figure 39. Competition for nutrient diffusing in at x = 5 cm. ........................................... 94 Figure 40. Competition between Pseudomonas KC and E. coli in swarm plates. ............. 97 Figure 41. Colony counts in competition experiment. ....................................................... 98 ix Figure 43- 34' 652111643- 5* Figure 44- 5- figure 45. Si Figure 46. .\1 figure 47. .\I Figure 48. CI Figure 49. Ni Figure 50. .\': Figure 51. P: Figure 53. P: Figure 53.. .\-1 Figure 54. C Figure 55. Sr Figure 56. Sr Figure 57. P: Figure 58. C Figure 59. T Figure 42. Mutant selection by chemotactic response (source at top of pictures). .......... 100 Figure 43. E. coli responding to chemoattractant gradient. ............................................. 104 Figure 44. E. coli responding to inhibitor gradient and chemoattractant gradient. .......... 107 Figure 45 . Simulation of experiment shown in Figure 42. .............................................. 109 Figure 46. Mutant appears far back in population. .......................................................... 111 Figure 47. Mutant now chemotactic to S. ....................................................................... .112 Figure 48. Chemotactic rings of Pseudomonas KC at varying nitrate levels. ................. 116 Figure 49. Nitrate concentration = 50 mg/l, varying acetate ............................................ 117 Figure 50. Nitrate concentration = 500 mg/l, varying acetate. ......................................... 118 Figure 51. Pseudomonas KC responding to glucose (left) and aspartate(right). ............. 119 Figure 52. Pseudomonas KC responding to acetate and nitrate in a DGC. ..................... 122 Figure 53. Movement of wave around obstacles in swarm plates. .................................. 123 Figure 54. Chemotaxis of Pseudomonas KC around an object. ...................................... 125 Figure 55. Set-up to make glass—bead wall in DGC. ........................................................ 127 Figure 56. Swarm plates with various chemoattractants .................................................. 132 Figure 57. Pattern formation in the DGC ......................................................................... 133 Figure 58. Chemotaxis through a hollow object. ............................................................. 134 Figure 59. The alternating-direction implicit ethod. ........................................................ 135 x nutrie chem chem diffus diffu.t dIIfUr arera glucc IIULFIT cell 1 cell 1 cell 1 inhit mass mass flux flux fluX total inhil LIST OF NOMENCLATURE nutrient saturation constant chemoattractant Q saturation constant chemoattractant S saturation constant diffusion coefficient of H diffusion coefficient of Q diffusion coefficient of S average maximum flux glucose concentration nutrient density cell flux cell receptor saturation constant for chemoattractant Q cell receptor saturation constant for chemoattractant S inhibition constant mass transfer coefficient for Q mass transfer coefficient for S flux of cells due to chemotaxis to S flux of cells due to chemotaxis to Q flux of cells due to random motility total number of receptors for a chemoattractant inhibitor P density xi chem nucrc chern thne ceH d ' . rolur niaxi maxi speci chen perc yield cent ungl chen Chen rand rand Cher dyn: ena chemoattractant Q density microbiosensor reading chemoattractant S density time cell density volume of DGC arena maximum growth rate on H maximum specific growth rate modified for inhibition specific chemoattractant consumption coefficient for Q specific chemoattractant consumption coefficient for S chemotactic velocity percent oxygen saturation yield coefficient for cells growing on H cell tumbling frequency single cell swimming speed chemotactic sensitivity coefficient chemotactic sensitivity coefficient modified for inhibition random motility coefficient random motility coefficient modified for inhibition chemotactic response factor for chemoattractant S dynamic competition factor xii N ' goflossible for chemotaxis to be of a beneficial nature, and then to deveIOp ways to engineer hose systems so as to realize the full benefits offered by the chemotactic response. The hree applications that were studied were microbial competition, selection of mutants, and n situ bioremediation. Section 2 explores the two main objectives in more detail, and describes their ignificance to the scientific community. Section 3 gives background on the experimental .nd modeling methods employed by others to study chemotaxis. The results of this research 12101e methods thatr 57 describe I and the work srrnrnan 01 1h Include hell in the ma Trier are incl chemotactic 11 outlines the alt equations. Apr Which is prese program. The I the quantities p the model is de 2 search project are presented in the next four sections. Section 4 examines the tools and :thods that were developed to study chemotaxis over the course of the project. Sections 1 describe the applications that were identified for improvement through chemotaxis, i the work that was done to engineer the systems to make use of chemotaxis. A nmary of the work, and the conclusions that were obtained are given in Section 8. Included in Appendix A are results of some interesting experiments that did not fit 11 in the main text of the thesis, or have not been fully analyzed for their importance. ey are included to show the range of growth patterns that can arise due to the :motactic response, including rings and geometric patterns. Appendix B briefly :lines the altemating-direction implicit (ADI) algorithm for solving partial differential Iations. Appendix C gives instructions for the FORTRAN ADI model, one version of ich is presented in Appendix D. Appendix E gives the input files for the FORTRAN gram. The Matlab files used to produce many of the graphs and to calculate some of quantities presented in the thesis are given in Appendix F. A dimensionless version of model is developed in Appendix G. 2. OBJEC' 3.1 DeveIOj The f1 modeling met The diffusion used in this re Emerson et at SIstem. Cons different oper responding 1. chemoattracra inhibitor grai quantitative 1r without actuaj T0 pn andmicroorg; Ihesfiltaramet the Iager diffr; OilherandOm OBJECTIVES AND SIGNIFICANCE . Development of tools The first objective of this research was to improve existing experimental and deling methods for studying chemotaxis, or to develop new methods, as necessary. e diffusion gradient chamber (Emerson et al., 1994) was the main experimental tool :d in this research. The DGC was well established and tested for chemotaxis studies by Ierson et al., but some further enhancements to the method were made to optimize the :tem. Considerable research was devoted to developing mathematical models of ferent Operating modes of the DGC system. These modes included one cell population .ponding to one chemoattractant; two cell populations responding to multiple :moattractants; and two cell populations responding to two chemoattractants and an Libitor gradient. These models were experimentally validated and then used for mtitative interpretation Of experimental results, and for exploring chemotactic behavior :hout actually running time-consuming experiments. TO produce useful simulations that accurately modeled real chemical gradients I microorganisms, parameter values had to be Obtained. Several methods of measuring 36 parameters were developed as part Of the research effort. These methods included laser diffraction capillary assay (Schmidt et al., 1997) which allowed for measurement he random motility coefficient, and microelectrodes and microbiosensors (Peteu et al., 5), Which enabled high—resolution measurements Of gradients of oxygen and glucose. 2.2 Applie The St chemotactic n has the first s of bacreria c competition 1 hipothesized position them The DGC pror methods to \‘a Anom- mutant strains to draw muta Concentration 11161th 0f SE obtained for [1 A thin The bacteriun FIOdchng 11 {if cllrrentty bein Applications for chemotaxis The second Objective Of the research was to identify applications in which the totactic response could be exploited to benefit the application. Microbial competition he first such application selected. Many studies have examined how two populations acteria compete in well-mixed environments, but much less is known about petition in spatially structured systems, where chemical gradients can occur. It is tthesized that chemotaxis can impart a competitive advantage by allowing cells to :ion themselves in conditions optimal for their survival. The mathematical model of )GC provided a means to explore this hypothesis. Preliminary work on experimental rods to validate the model's predictions was also initiated. Another application benefitted by the chemotactic response was the selection of [III strains of bacteria that exhibited a desired trait. Specifically, chemotaxis was used raw mutant strains of Escherichia coli into a region of the DGC having a higher entration of an inhibitor than their non-mutagenized progenitors could tolerate. This Od of selection by chemotaxis resulted in an improved strain of E. coli being ned for the industrial fermentation of a valuable product. A third application for which chemotaxis was studied was in situ bioremediation. )acterium Pseudomonas KC, which is able to degrade carbon tetrachloride without Icing harmful by-products, was chosen as the model organism. Pseudomonas KC is My being used in Schoolcraft, MI, to bioremediate a contaminated aquifer. The chemotactic 1 medium on m otactic behavior of Pseudomonas KC was studied, and the effects of a porous um on motility and transport were observed. 3. CHEMI 3.1 The ch hr mohhu' thnnnella C 193.1. “11611 he flageha r hhenthe ceI mmhmdn ofthn senes hgmel. Inanenviron the Presence IEMOTAXIS BACKGROUND he chemotactic mechanism otility and chemotactic response Of the enteric bacteria Escherichia coli and rella typhimurium have been studied in detail (Berg and Brown, 1972; Macnab, When these peritrichously flagellated bacteria spin their flagella counterclockwise, gella move in a synchronized bundle, causing the cell to swim forward, or run. the cell reverses the spin of the flagella, the bundle separates, causing the cell to and randomly reorient itself before running in a new direction. The overall effect series of runs followed by tumbles is referred to as a random walk, as shown in 1. Figure 1. Random walk .vironment with no chemoattractants, the run length is independent of direction. In ence of a chemoattractant gradient, the cells monitor the time rate of change of occupancy t concentratior extended run frequency ap result is a bin in Figure 3. 3'2 EXperi EXper assay (Adler, 1972; NOSSal (SFDCI (Fort 7 upancy of receptors for the chemoattractant during a run. If an increasing centration gradient is encountered, the tumbling frequency decreases, causing an ended run length. If a decreasing concentration gradient is encountered, the tumbling uency apparently remains at the basal level (Berg, 1988; Macnab, 1987). The net tlt is a biased random walk toward the higher chemoattractant concentration, as shown igure 2. Low I I . I : High; Concentration Figure 2. Runs and tumbles in the presence of a gradient. Experimental systems for measuring chemotaxis Experimental systems used to study microbial chemotaxis include the capillary (Adler, 1972; Nikita er al., 1992; Staffeld et al., 1987), the motility plate (Adler, ; Nossal, 1972; Wolfe and Berg, 1989), and the Stopped Flow Diffusion Chamber C) (Ford, 1992; Ford et al., 1990). In the capillary assay, the open end of a capillary filled with bt for a fixed at entered the 1 thiferent co mathematical Iauffenburge random motil A mo hhing a unif of the plate 1 gradient is To The motility Cellular met; Chtrnotactic motility plate TheS are introducer 01 interegr u appIOXImated an initially Sit by Iight mun :d with buffer solution containing a chemoattractant is inserted into a cell suspension a fixed amount of time. The capillary is then removed, and the number of cells that _.> fired the capillary are counted. By comparing this number for capillaries having Lrent concentrations of chemoattractant, chemotaxis can be quantified. A ematical model of the capillary assay has been developed (Rivero-Hudec and enburger, 1986) that allows calculation of two important modeling parameters, the om motility coefficient and the chemotactic sensitivity coefficient. A motility plate consists of a petri dish containing a semi-solid agar medium ng a uniform distribution of a consumable chemoattractant (Adler, 1972). The center le plate is inoculated with cells. As they grow and consume the chemoattractant, a ient is formed that induces migration of the cells outward from the inoculation zone. motility plate is a simple system to use, but, because the gradient is formed by tlar metabolism, the gradient is difficult to control and quantify. Also, the iotactic response to non—metabolizable chemoattractants cannot be studied in a city plate. The SFDC is a rectangular chamber into which two impinging streams Of medium atroduced. One stream contains the chemoattractant, and the other contains the cells terest. When the flows are suddenly stopped, a step—change in chemoattractant is gximated where the two streams intersect. Diffusion of the chemoattractant creates (ft 11 <2 :itially steep gradient that decays over time. The chemotactic migration is measured h ght diffraction. A mathematical model of the system has been developed to calculate the chemoatt cells. .\'one tell characte hate develop microbial che previOUSIy PU 9 chemoattractant gradients and to analyze the resulting chemotactic mi gration Of the s. 1 None of the experimental systems described above are convenient for establishing 1 characterized, steady-state gradients or multiple gradients in multiple directions. We 3 developed the Diffusion Gradient Chamber (DGC), shown in Figure 3, for studying lrobial chemotaxis under such conditions. Details of the DGC system have been Reservoir j n 11 l n, e _________________ I e 4 ET Arena ("m a 2 er— ‘5: ...... E _ __'. _ gm ES - . _________________ , ’t e 2 t 9 Gasket E1 1 U El Semi— ermeoble / \ mem rane Inlet Outlet Figure 3. Diagram of the diffusion gradient chamber. ously published (Emerson et al., 1994), so only a brief description is given here. The consists of a square arena (5 cm x 5 cm x 1 cm) bounded by a reservoir on each «in? side. Each re byan imperr arena and res gel through chemoattract concentratior concentratior across the gr then the char monitored lit The 1 31111111316 gm Chamber. GT inoculation, Simultaneous Tilt DGC ha microsengOr inOCUIEITIOn i“Oculation i maintain an . gaseous reag 10 side. Each reservoir is separated from the arena either by a semi-permeable membrane, or by an impermeable silicone sheet, depending on whether substrate diffusion between the arena and reservoir is desired. The arena is filled with medium containing dilute agarose gel through which the cells can readily swim. Different concentrations of the chemoattractant(s) are maintained in the reservoirs. Diffusion from the higher concentration reservoir (the source) through the membrane and across the gel to the lower concentration reservoir (the sink) results in the development of a continuous gradient across the gel. The gradient is allowed to establish for a specified amount of time, and :hen the chamber is inoculated. Growth and movement of the microbial populations are nonitored from above by light diffraction. The DGC has several advantages over other methods of studying chemotaxis. vIultiple gradients in different directions may be established simultaneously in the :hamber. Gradients may be allowed to approach a linear steady-state profile before noculation, or transient gradients can be used. These gradients can be initiated imultaneously, or in a staggered fashion to simulate pulsed influxes of chemoattractants. The DGC has been designed with a removable lid so that samples may be withdrawn or microsensor readings taken without sacrificing the experiment. Many different moculation protocols may be used, including uniform inoculation across the gel, loculation in a line, or at a point. Gas ports in the chamber provide the possibility to laintain an inert gas headspace above the arena for anaerobic experiments, or to supply aseous reagents to the microbes. 4. DEVEI 4.1 lmpro 4.1.1 Image In prr the DGC wet was laboriot course of [In were detelo; The f rcx ccnc The Pulnix c ‘Qllanta, Mc Written by Se Chemotactic Smdiertarbeit PICiures at pp 4‘0 for Wind Two ca[Horas (C01 other to a P 0‘ DEVELOPMENT OF TOOLS Improvements to the DGC method 1 Image acquisition system In previous work done in the DGC system (Emerson et al., 1994), photographs of DGC were taken manually with a 35 mm camera. This method gave good results, but laborious and did not lend itself well to computerized image analysis. Over the .‘se of the work presented in this thesis, several automated image analysis systems :developed. The first computer-based image analysis system was built using a PULNiX TM- 1 CCD-camera (PULNiX America Inc., Sunnyvale, CA) mounted above the DGC. Pulnix camera was attached to a PC through a WinVisionPro video capture board tnta, Mountain View, CA). A time-lapse capture program called AutOCap was en by Sebastian Schmidt ("Development of Novel Methods to Measure Random and notactic Microbial Motility at the Community Level", Sebastian Schmidt, ienarbeit, Michigan State University, 1995) and used to automatically capture res at pre-set intervals. Captured images were converted into RAW format in Matlab )r Windows, and then analyzed. Two other image analysis systems were installed based on Color QuickCam ras (Connectix, San Mateo, CA). One Quickcam was connected to a PC, and the to a Power Macintosh 7200/90. The software included with the Quickcams allowed 11 for timed-in PL'LViX sys afreetrare ir 4.1.2 Quan To q curve relatir. has grown glycerol for were taken madmac acquired. Tl (lithe cultur bar. Anothe medium, 0p added [0 1111 IHOCulated, ' The 1116 1:100,()( average but 12 ned-image capture. Images from these units, and, later in the research, from the liX system, were analyzed using NIH Image (http://rsb.info.nih.gov/NIH—I1VIAGE/), ware image analysis program. Quantification of cells by grayscale analysis To quantify the number of cells in a DGC through image analysis, a calibration relating cell number to image grayscale was developed. Escherichia coli HCB 33 ;rown in M63 medium (described in more detail in Section 4.2.3) with 5 mM 01 for 30 hours at 30°C. Replicate measures of the optical density of the culture taken at a wavelength of 630 nm. Forty milliliters of the cell culture were then to a clean DGC containing no agar, and three consecutive images of the DGC were ed. The DGC was then emptied and cleaned for the next reading. A 0.1 ml aliquot culture was placed in triplicate on LB agar plates and spread uniformly with a glass nother volume of the cell culture was then diluted by a factor of 1:1.5 with M63 m, optical density measurements were taken, 40 ml of the diluted solution was to the DGC, three more images were acquired, and three more LB plates were Ited. This process was repeated seven times at the dilutions shown in Table 1. The plates from each dilution were incubated at 30°C for 18 hours. At this time, 00,000 dilution plates had the most suitable number of colonies to count. The a number Of colonies in each plate was 448, with a standard deviation of 36 s. The original culture therefore contained approximately 4.48><108 cells/ml. Table 1. Dilution factor and OD for grayscale calibration 13 Dilution Optical Density (ml cells/ml total) (7t=630 nm) no dilution 0.762 0.789 1:1.5 0.559 0.577 1:2 0.453 0.424 1:3 0.327 0.299 1:4 0.322 0.270 1:5 0.191 0.192 1:10 0.082 0.066 1:15 0.054 0.041 1:20 0.025 0.028 1:100 -- -- 1:1000 -- -- 1:10,000 -- -- 1:100,000 -- -- vlatlab 4.0 for Windows was used to analyze the images. The image was :d to RAW format with Image Alchemy, and then read into Matlab as an image mefleme imgfima dilution. Fin: tononnahze meraged at ihmmdm Tm DGC.Thn wfihiyhr. [his (IISSen mmmmw 14 . The average grayscale value of each picture was calculated in Matlab. The three .ges for each dilution were then averaged together to give one grand average value per ttion. Finally, each grand average was divided by the maximum of the grand averages tormalize the values into a range from 0 to 1. The Optical density readings were also raged at each dilution. The relationship between optical density and grayscale is strated in Figure 4. Normalized grayscale 0 0.1 0.2 0.3 0.4 0.5 0.6 Optical Density Figure 4. Normalized grayscale vs. Optical density for E. coli. The normalized grayscale was next plotted against the number of cells in the C. This yielded a linear calibration, as shown in Figure 5. Other work in collaboration 1 Tyler Ames, a Ph.D. candidate in Dr. Worden's research group, but not included in dissertation, has also shown that soybean plant cell concentration has a linear tionship with grayscale in the DGC (Ames, 1997). 4.1.3 Metht An e the membra following Se arena by a s Was iTXtd in [01116 grour harden, The take Place in ahOFiZOrltal A su 3000 “1L of 15 1 9 0.8 ~- 8 9 <6 -_ a 0.6 .0 a g 0.4 ~- E o 2 0.2 0 l l l 7 O 1 2 3 4 5 # Of cells x 10'8 Figure 5. Calibration curve for E. coli in DGC. 3 Method to measure membrane mass transfer coefficient An experiment was designed to measure the mass transfer coefficient, ks, across membrane. This parameter is important in the modeling work introduced in the wing sections. A DGC system was prepared having one reservoir separated from the a by a semi-permeable membrane, and the other three sealed. After the membrane fixed in the chamber, the chamber was stood on one side with the membrane parallel 1e ground. 400 11L Of 0.3% agar were spread onto the membrane and allowed to en. The agar layer was applied to duplicate any fouling of the membrane that could place in a DGC experiment. After the agar layer solidified, the DGC was returned to rizontal position, and filled with 45 ml of RO water. A sucrose solution was pumped through the reservoir from a flask containing lmL Of 30 g/L sucrose. The solution was recycled back into the 3000 mL flask. A stir bar was plac model for thi where S is Il A is the are. concentratio 5:0 and line Sugar cone. measuremer 7’ _‘ ix.) In (Ssource/(Ssourco-S)) O _(D 51>. or) 0% 16 ;ar was placed in the DGC, and the DGC sat on a magnetic stir plate. The mathematical odel for this system is l V flicks/us —S) (1) arena source Where S is the concentration of sugar in the DGC arena, Varena is the volume of the arena, x is the area of the membrane (3.06 cmz) available for mass transfer, and Ssource is the oncentration in the source flask. Equation ( 1) can be solved with the initial condition =0 and linearized as —S V source arena 2 1n[ Ssource J: kSA t ( ) S . ugar concentrations were measured over time by HPLC analysis. The results of leasurements for both sucrose and glucose are shown in Figure 6. Sucrose Glucose 1.6 1.6 53 1‘2 k=o.28 chh :3 1-2 k=0.31 cm/h 25 08 - £3: 03 l E 3 R2 = 0.9988 3 R2 = 0.9977 ‘5? 0-4i 9’; 0.4 l .5 2 0 I l —*r 0 cr 1; l l O 20 4O 6O 80 O 20 4O 60 80 time(hr) time(hr) Figure 6. Determination of the membrane mass transfer coefficient. 4.1.4 Pastel In St conventiona' arena. \l'hel process beg mathematicz state diffusit Where 3 is 1 Equation l 3 Where \ 13 boundary CO 17 f4 Faster approach to steady-state gradients ! . In some instances, steady-state gradients were desirable in the DGC. In the Eventional Operating mode, gel without the gradient chemical was poured into the na. When the gradient chemical was introduced into one reservoir, the diffusion Fess began, and eventually a linear steady-state gradient would be established. A hematical model of the system was developed to describe this process. The unsteady- e diffusion equation, for the case with no cells present, is given by 35 ——- = D V25 8: 5 (3) :re S is the gradient chemical and D3 is the diffusion coefficient. In two-dimensions, .ation ( 3) can be rewritten as 825 825 (4) a = . 55*? ;re y is the direction parallel to the gradient and z is the vertical direction. The dary conditions at the membranes are given by I ( 5) ’35 = 1(5— (SIFO " SSink) ay D. he sink reservoir (position y=0 cm) and ( 6) (3‘ —k $16 = —D-: (S|y=5 — SSource) hhwm assumed to boundary CO The equation iii eoui'entiona steady-state. state solutic solution at 1 within 95‘} necessary tc A ne achieve the Work refere schematical COttcentratic Wired and “0 gradient 50% 0f the 18 ' the source reservoir (position y=5 cm). The source and sink flask concentrations were ;umed to be constant throughout the experiment. For all other boundaries, the no-flux undary condition was applied l sts = o (7) The Matlab Partial Differential Equation Toolbox was used to solve the diffusion 1ation with a finite element method. The predicted time to reach steady—state for the aventional method was 33 days, using sucrose as the model chemical. To define ady-state, the elliptical equation was solved by the tool box to give the actual steady- te solution. Then the parabolic equation was solved, and compared to the elliptical ution at several time points. The first time point at which the parabolic solution was bin 95% of the elliptic solution at every node in the DGC was defined as the time :essary to reach steady-state. A new method of pouring gel slabs was developed to reduce the time necessary to .ieve the steady-state gradient. This method was developed as part of the plant cell rk referenced in Section 4.1.2. The new method of pouring the slabs is shown ematically in Figure 7. The DGC was tilted at an angle and Layer 1, containing a centration of the gradient chemical equal to that used in the source reservoir, was Ulmlll fitted and allowed to solidify. The DGC was returned to level, and Layer 2, containing h‘Illl “gradient chemical, was poured and allowed to solidify. Finally, Layer 3, containing 6 of the gradient chemical concentration used in Layer 1, was poured. The Matlab distance frorr 19 _U‘r—— 50m —————>l 0.4 cm ayer nocu um T T Layer 2 2.2 cm 2 Layer 1 L, Figure 7. Three-slab pouring method. Matlab simulation of this system, depicting the concentration as a function of tnce from the sink reservoir, is shown in Figure 8 for three time-points. '20 10 0 Sucrose (g/l) 5 hr 012345 Position (cm) Figure 8. Simulation of approach to steady-state. :ime necessary to reach steady—state for the three—slab pouring method is predicted to 1 hours, which is approximately 20 times less than the conventional pouring method. 4.15 Recy Ano flasks. In t through the method wet the flasks a flask m concentratit lift/week result is spt {185k conce mode. “5 Micr Mic. Worden 10 1 01.. 1996a; microelecm "HCFOdectn conslimes g miceoelecm 20 5 Recycle of source and sink flasks Another in the operating method was the use of recycle in the source and sink ks. In the original DGC method, fluid is pumped out of the source or sink flask, >ugh the correSponding reservoir on the DGC, and then to a waste container. This thod worked well, but the time an experiment could run was limited by the volume of flasks and the pump speed. The idea of recycling medium from the reservoir to the ;k was tested. Sucrose was the gradient chemical, at an initial source flask reentration of 20 g/L. In six different experiments, only 0.26i0.13 g/L-week (or %/week) were lost from the flasks due to transfer across the membrane. Although this ult is specific to the plant cell system and to sucrose, it suggests that in general, the 3k concentrations will not change significantly if the system is operated in recycle de. .6 Microsensors and microbiosensors Microsensors were developed by David Emerson, Serban Peteu, and Mark rden to measure oxygen, glucose and other chemicals (Peteu et al., 1996; Emerson et 1996a; Emerson et al., 1996b). Briefly, the microsensors are Clark-type oxygen roelectrodes. The addition of an enzyme, such as glucose oxidase, to the tip of the oelectrode produces a microbiosensor. The enzyme catalyzes a reaction that umes glucose and oxygen. The rate of disappearance of oxygen is measured by the >electrode, and allows for the glucose concentration to be measured. A P< for the glue linked to th hmmm function of lit the data I Where R is percent oxy least-square Parameters calibration l 21 A portion of this research project was devoted to deve10ping a calibration method the glucose microbiosensors. Because the reaction at the tip of the microbiosensor is ked to'the oxygen concentration, calibrations at varying oxygen concentrations had to performed. A group of microbiosensor readings, in picoamps, was obtained as a action of percent oxygen saturation and glucose concentration. The simplest model to the data to is a linear plane, whose equation is given by R=aG+bX+d (8) here R is the microbiosensor reading, G is the glucose concentration, and X is the :rcent oxygen saturation. The Solver tool in Microsoft Excel 7.0 was used to find the ast-squares best fit between the experimental points and the model, by varying the trameters a, b, and d simultaneously. Four isoclines of the resulting three-dimensional libration plane are shown in Figure 9. An microbiose; dmuhaneOI medlUm W; for m0f€ dE bacteria hat 22 400 00 mM glucose 0 A136 mM glucose 2 __ x476 mM glucose 9; 300 [1816 mMglucose 8’ 1.; Symbol: experiment 0 - 9 200 “c Line = best fit plane X 8 8 3 ‘ II a 100 —— . A " ° A El 0 l + l #1 r O 20 40 60 80 100 % 02 saturation Figure 9. Three dimensional glucose microsensor calibration plane. An example of an experiment in which the oxygen microsensor, the glucose icrobiosensor, and the image analysis system for cell grayscale were used nultaneously is shown in Figure 11. A swarm plate containing 1 mM glucose in M63 :dium was inoculated with 10 {AL actively growing E. coli HCB 33 (see Section 4.2.3 f more details on medium and strain). After approximately 24 hours of growth, a ring of :teria had formed, as shown in Figure 10. klicrosensc An image V 11 Shows chemotactit 23 Figure 10. E. coli forming ring in response to glucose. Microsensor and microbiosensor readings were taken across the ring at a depth of 1 mm. An image was recorded of the swarm plate with the Pulnix image analysis system. Figure 11 shows the combined results of the sensor readings and the image analysis. The chemotactic wave (see Section 4.2) is visible where the gradients have the steepest slopes. Figure l 4.2 Den The Coupled CC or Carbon 3 Where Ju i lllllfiem (1‘ review, St migration 24 1200 30 _ c‘a" ”70 100° ° n “ ,. Edge of wave A t‘ _ 60 % O E 800 -— .. ii a "’ - 50 9 _o' .. <5 C H o .0 8 600 " __ 40 g ‘1’ N 8 n __ 30 O .3 400 .- +Glucose g ‘9 +oxygen .. ' Inoculation__ 20 o\.. 200 __ —Cell Profile 90"“ - 10 U 0 . . “_-..___: : : :: : : , , 0 -15 -1o -5 0 5 10 15 20 Distance (mm) Figure 11. Use of microsensors, microbiosensors, and image analysis, DGC modeling 4.2 Development of model The mathematical model of the Diffusion Gradient Chamber system consists of coupled conservation equations for the microbes, the chemoattractant(s), and the nutrient or carbon source. The cell balance equation may be written as a (9) §=—V-Ju +f(H)u where Ju is the cell flux, 11 is the cell density, and f(H) is a function for cell growth on a nutrient (H). Many forms of the flux equations for chemotaxis have been suggested (for a review, see Ford, 1992). A commonly used constitutive equation to describe cell migration at the population level is that derived by Keller and Segel (1971): there it is model pres chemotactir there \"',,5 chemotactir The the RTBL r W has be IItodel for assumes th; to the Spa (Flymjer, e. 25 J, =—uVu+Vuu (10) where u is the random motility coefficient and Vu is the chemotactic velocity. For the model presented here, the Keller-Segel equation has been modified by adding a second chemotactic flux, to give Ju =—,uVu+VuSu+V,Qu (11) l where Vus is the chemotactic velocity in response to a chemoattractant S, and VUQ is the chemotactic velocity to a chemoattractant Q. Combining Equations (9) and ( l 1) gives au 8: (12) =uV2u—V-(V.su)—V-(V.Qu)+fu The chemotactic velocity equation prOposed by Rivero, et al. (1989) (known as the RTBL model) is strictly applicable only to one-dimensional movement of the bacteria, but has been shown to give good agreement with a more rigorous three-dimensional .nodel for a range of parameter values (Frymier, et al., 1992). The RTBL model also assumes that the temporal gradient of chemoattractant has a small contribution compared :0 the spatial gradient, and evidence to support this assumption has been provided :Frymier, et al., 1994). The RTBL chemotaxis term is given by (13) VS =vtanh w—MVS} “ (K0, +5) where t) is total numbr for the rec velocity car where 105. terms of llll A similar e Where X00 dissoeiatio Chemotacti from man that both c assay 26 here 1) is the swimming speed of a single cell, 0 is the tumbling frequency, NTS is the ital number of receptors for the chemoattractant S, and KDS is the dissociation constant tr the receptor-S complex. For shallow gradients, the modified RTBL chemotactic :locity can be approximated by KDS VS (14) (Km + Sf Vus = 105 here X08: the chemotactic sensitivity coefficient to the attractant S, can be expressed in arms of individual cell parameters by X05 =002Nrs (15) similar equation can be written for the chemoattractant Q as K (16) DQ 2 VQ KDQ+Q) VuQ = log ( here XOQ is the chemotactic sensitivity coefficient to the attractant Q and KDQ is the _ ssociation constant for the receptor-Q complex. Rivero et al. (1989) showed that the remotactic sensitivity coefficient and the random motility coefficient can be calculated om individual cell parameters, and Rivero-Hudec and Lauffenburger (1986) showed at both can be measured experimentally through population assays such as the capillary say. The and (16) w yield 6+“er dt where v is half-saturat chemotactit better thar reproducet achieved g coli to simt The Monod—[yp Chemoattra “Puke. Th nutrient (l chemoattra The 27 The cell balance used for the DGC system is derived by combining Equations (14) d (16) with the random motility flux term and a term for growth (on substrate H) to :1d 9a K K vH (17) _= V2 _ 0 —L— _ . DQ 9t ’1 u ZOSV [[(Kos + S)2 }VS] XOQV [[(KDQ +Q)2 }VQ]+ C0 + H u tere v is the maximum specific growth rate of the cells growing on H, and C0 is the lf-saturation constant for growth on H. The assumption of simple additivity of the :motactic responses to multiple attractant gradients has been shown to be as good or :ter than more complex interaction models, but lacks the ability to completely iI‘OdllCC experimental observations (Strauss, et al., 1995). Boon and Herpigny (1986) rieved good agreement with experimental results when modeling the response of E. 7i to simultaneous gradients of glucose and oxygen using this assumption. The cell growth term and the chemoattractant consumption terms are modeled as mod-type (Bailey and Ollis, 1986) saturation processes. Growth of cells due to :moattractant uptake is assumed to be negligible compared to growth due to nutrient .ake. This assumption is reasonable for many experimental systems because the rient (usually glycerol) is present at a much higher concentration than the mo attractant . The two chemoattractant balances are given by 1255 u (18) as 2 ———= VS- 31‘ D5 CS+S and where D; a chemoattra consumptit The Where DH ; 0n H. In 1 diffusion c aniougly 4-2-1 Bou A 1' illustrate tl 2'5 cm, am on all sou, 28 vQQ (19> where D5 and DQ are the chemoattractant diffusion coefficients; vs and VQ are the specific chemoattractant consumption coefficients; and Cs and CQ are the saturation constants for consumption of S and Q, respectively. The nutrient balance is given by a—H-=DHV2H— 3t vH L (20) C0 + H YH where DH is the nutrient diffusion coefficient, and YB is the yield coefficient for growth on H. In both of the chemoattractant balances and in the nutrient balance, the respective diffusion coefficients are assumed to be constant and the medium isotropic. We have previously confirmed that substrate diffusion into the DGC is accurately described by the model (Emerson 81‘ al., 1994). 1.2.1 Boundary conditions A two-dimensional schematic representation of the DGC is shown in Figure 12 to ,llustrate the geometric parameters used in the mathematical model. The dimension R is 2.5 cm, and r is 2 cm. For the cell balance, a zero total flux boundary condition is applied in all boundaries (Q) of the chamber. This boundary condition, Shown below, states that the overall [02610 on 1 Net boundaries diffuse ac transfer. T] where Sm resfirt'oirs, Penneable Ou Walls of [h confllilOn 29 3 overall flux, given by the sum of the random motility and chemotaxis fluxes, is equal zero on the boundaries. 2 KDS KDQ _ (21) {#V u ‘Zosv'[[(1{m +S)2 lVSl—ZOQV-lllKDQ +Q)2 }VQlln - O Neumann boundary conditions are used for the chemoattractant S on the undaries where the reservoir is open to the arena. The chemoattractant is assumed to fuse across the semi—permeable membrane, which imposes a resistance to mass nsfer. This boundary condition is written as ( 22) .ere Sres,S and Stem are the concentrations of the chemoattractant in the south and north ervoirs, respectively, and ks is the mass transfer coefficient for S across the semi- meable membrane. Outside the reservoir openings, there is no flux of chemoattractant across the 118 of the chamber: 35 (23) 85 __ :0 832% 3y =0 and 1‘0 :he reservoir is sealed with a non—permeable membrane, then the no flux boundary tdition would apply across the entire side. The boundary conditions for Q and H are similar ma thesis. the initially pr: C 01le Ilil at. 4.2.2 Init h inoculated shape of ti Where up, i: mdw'we exPerimen lhal X':0 3 Th initial cop Concentrat XWME( 413E“ Ex Sirepmmy‘ 3O irnilar mathematically to those of S. For all simulations presented in this section of the hesis, the concentration of Q in each Open reservoir was equal to the concentration of Q nitially present in the arena, and the concentration of H in each reservoir was equal to the oncentration of H initially in the arena. ..2.2 Initial conditions In the experiments used to validate the model, the center of the chamber was aoculated with a micropipette. An exponential function was used to approximate the hape of the injected cell peak: u (24) 11(X’, I): 0 y exp(w1/x’2 + y’z) rhere uo is the initial concentration of cells and w is a peak width factor. The values of uo nd w were calculated to yield the same number of cells in the peak as were added xperimentally (~4x107 cells, Emerson et al., 1994). The variables x' and y' are defined so rat x'=0 and y'=0 occur at the center of the arena. The initial condition for the chemoattractant S is S(x,y)=0 for all x and y. The titial condition for H is that H(x,y)=Ho for all x and y, where H0 is the initial )ncentration of H in the arena. The initial condition on Q is that Q(x,y)=Q0 for all x and . where Q0 is the initial concentration of Q in the arena. 2.3 Experiments used to validate model Experiments in the DGC were carried out as described in Emerson et al. (1994). A reptomycin resistant strain of E. coli HCB 33 (=RP437) that is Wild tYP€ for motility and chemoraxis w This medium was methionine. and t] Tne glycerol ser chemoattractant i and either 0.1 m! 1163 medium on initially containe. in LB broth. In thesee north the sink. 1 pm pore-size p0] 31 chemotaxis was used. A mineral salts medium (M63) was used for all experiments. 3 medium was supplemented with 2 mM glycerol, the amino acids histidine, leucine, hionine, and threonine required for this auxotrophic strain, and 125 ug/l streptomycin. 2 glycerol served as the carbon and energy source for growth, but it is not a moattractant for E. coli. The source reservoir contained supplemented M63 medium either 0.1 mM or 1 mM aspartate, while the sink reservoir contained supplemented 3 medium only. The medium in the arena was stabilized with 0.15% agarose and ially contained no aspartate. The E. coli were grown overnight with aeration at 30°C .B broth. In these experiments, the south reservoir (see Figure 12) was the source, and the th the sink. These source and sink reservoirs were separated from the arena by 0.05 pore-size polycarbonate filter membranes (Poretics Corp., Liverrnore, CA). v Figure The cast silastic sheeting autoclaved at 12 in the flasks as I sterile filters. Th through a sterile Both the reservoirs at 2.5 Operated for a 32 North Reservoir '= rn O a 8 a, H a) 30 a: (D 01 H—R———>I g ‘8 ' 2 ' O E . =' I d v “‘f r . X South Reservoir Figure 12. Schematic diagram of DGC, showing modeling coordinates. The cast and west reservoirs were sealed off from the arena with non-permeable silastic sheeting (Dow Corning, Inc., Midland, MI). The source and sink liquids were autoclaved at 121°C to sterilize, and allowed to cool. To prevent a vacuum from forming in the flasks as liquid was pumped out, the flasks were vented to the atmOSphere through sterile filters. The gas port to the DGC (see Figure 3) was also vented to the atmosphere through a sterile filter, to allow the headspace to be replenished with fresh air. Both the source and sink solutions were pumped through their respective reservoirs at 2.5 nil/hr. To allow the aspartate gradient to partially form, the system was Operated for a given time before inoculation. The amount of time the gradient was allowed 10 form t the arena of the D the cells et'enl." [l 4.2.4 Computer The math: algorithm (Carnal difference equati equation. The fit only in the y-dire ADI method is a discretization err More details on t FORTRAN 77 at Q. and H were sr Was added, in or from the FORTE MATLAB progr: 4.2-5 MOdeling In order determined inde laser diffraction 33 allowed to form before inoculation was varied between experiments. The center point in the arena of the DGC was inoculated with 10 u] of E. coli using a micropipette to disperse the cells evenly throughout the depth of the agarose. 4.2.4 Computer simulations The mathematical model was solved with an Alternating Direction Implicit (ADI) algorithm (Carnahan et al., 1969; Chapra and Canale, 1988). The ADI method uses two difference equations to solve each two-dimensional unsteady—state partial differential equation. The first difference equation is implicit only in the x—direction and the second only in the y-direction. The equations are solved in succession at time steps of At/2. The ADI method is an unconditionally stable method with which convergence occurs with a discretization error of the order [(At)2+(Ax)2]. For the model presented here, Ax=Ay. More details on the ADI method can be found in Appendix B. The program is written in FORTRAN 77 and executed in the UNIX operating system. The balance equations for S, Q, and H were solved by the model for a specified amount of time before the cell balance was added, in order to simulate the gradient initiation time before inoculation. Output from the FORTRAN program was imported to MATLAB version 4.0 for Windows. The MATLAB program was used for image analysis and for graphical output. 4.2.5 Modeling parameters In order to validate the mathematical model, all of the parameter values were determined independently. The random motility coefficient, u, was measured using a laser diffraction capillary assay (LDCA) deve10ped in our laboratory (Schmidt et al., 1997). The LDC. movement in car and modeling. th indicate that the range of 0.15 to the experimenter. expenments. The assumed to be ir rate. in. was de battenal pOPUl'di? recalculate the g The diift correlation (Will and glycerol. re: 5111 by a eorrel (ire en. 1984). A “‘1 . tens growing on Values 1' dissociation "a ‘ . «Mil 11011: ‘I T - til‘ i ' pk liICIdl' 34 1997). The LDCA consists of a capillary tube through which a laser is directed. Cell movement in varying agar concentrations can be observed, and through image analysis and modeling, the random motility coefficient can be determined. LDCA experiments indicate that the random motility decreases linearly with agar concentration over the range of 0.15 to 0.30% agar. The random motility coefficient usedin the simulations is the experimentally measured value for the agar concentration used in the DGC experiments. The LDCA is described in more detail in Section 4.5. For this model, it is assumed to be independent of the chemoattractant gradient. The maximum cell growth rate, vH, was determined from the experimentally measured doubling time, rd, of the bacterial population using the relationship -1132 (25> v H t .1 to Calculate the growth rate. The diffusion coefficient for aspartate (D3) was calculated by‘the Wilke-Chang correlation (Wilke and Chang, 1955). The diffusion coefficients for Do and DH (oxygen and glycerol, respectively) at 25°C were obtained from. the literature, and adjusted to 30°C by a correlation for the temperature dependence of diffusion coefficients (Perry and Green, 1984). An approximate value for Ya was calculated by an electron balance for cells growing on glycerot. Values for the chemotactic sensitivity coefficient for aspartate, x05, and the dissociation constant for; the receptor-attractant complex for aspartate, K03, were taken Irom the iiterature. ‘s’alues tor the chemotactic sensitivity to various chemoattractants have been PUbli published by M3‘ The value range of those Ct specific chemoat the maximum ce csdmated by adj DGC iEmerson. glycerol). 8 (as; is given by kr-DC 0f the membrane and Chang. 195 predicted outcor coefficient. Valu 35 e been published (Ford, 1992). Values for the dissociation constant have been tlished by Macnab (1987). The values forthe saturation constants CH, Cs, and CQ were chosen to be in the ge of those commonly reported for E. coli (Bailey and Ollis, 1986). The maximum cific chemoattractant consumption coefficient (vs) was taken to be close in value to maximum cell growth rate on H. The mass transfer coefficients kH, kg, and kQ were mated by adjusting the k value for glucose, determined in previous work with the -C (Emerson, et al., 1994) by the ratio of the molecular weight of glucose to H ICCI'OI), S (aspartate), or Q (oxygen). The mass transfer coefficient across a membrane ;iven by k=Deff/l, where Deff is the effective diffusion coefficient, and l is the thickness he membrane (Cussler, 1984). Deg varies with the molar volume, V, to the -0.5 (Wilke , Chang, 1955) to ~06 power (Sitaraman et al., 1963). In results not shown, the :licted outcomes were relatively insensitive to small variations in the mass transfer fficient. Values for parameters used in the modeling simulations are given in Table 2. 416 Results 0 The cell those Observed “’38 introduced chamber at 21 cc Center of the ch. Thfi ext Figllre 13 The 36 Table 2. Parameter values for validation of model Parameter Value ll 0.0010 cm2 h" m 0.080 cm2 h" m 0.085 cm2 11'1 C0 4.08x10'6 gH cm? CQ 6.70x10’8 gQ cm'3 Cs 5.50x10'rgscmj’r DH 0.01 cm2 11'1 DQ 0.09 cm2 h"1 Ds 0.033 cm2 h'1 KDQ 3.30x10‘5 gQ cm' KDS 2.00x10'6 gs cm‘3 v 0.35 h" VQ 0.02 gQ gu"l h'1 vS 0.60 gs g," h“ YH 0.50 gu/gH 6 Results of validation simulations The cell growth and migration patterns predicted by the model were compared to ;se observed experimentally under two sets of conditions. In Run 1, 1.0 mM aspartate ,- MM 3 introduced into the south reservoir. Glycerol was initially present throughout the ) ‘lll‘ Timber at a concentration of 5 mM. After the gradient was established for 6 hours, the iter of the chamber was inoculated with E. coli. The experimental photos and modeling results at four time points are shown in Lure 13. The results of the computer simulation are shown in the left-hand column and photographs of tl modeling simulal observed experirr rate of cells. chz developed and in increased with ti approached the t source slowed 11 source. The thii chemoattractant In Run lOlHLVI aspanat in the left-hand exPeriment, the RU“ 1. chemota. “'3th were evit shapes of the 1 nglb pattern aPproached the 37 tographs of the experiment in the right-hand column. In both the photos and the fdeling simulations, lighter areas correspond to higher cell density. Several trends :erved experimentally were reproduced by the model. The first trend was the band or ire of cells, characterized by a locally high cell concentration near the pattern edge, that :eloped and migrated toward the chemoattractant source. The height of the wave peak [eased with time. The second trend was that the velocity of the wave decreased as it roached the chemoattractant source. Because the edge of the pattern closest to the rec slowed more than the rest, the pattern tended to broaden as it approached the rce. The third trend was i the less prominent wave that migrated away from the moattractant source reservoir. In Run 2, the chemoattractant concentration in the source (south) reservoir lmM aSpartate) was one tenth that in Run 1. Figure 14 shows the computer simulation ;he left-hand column and experimental photographs in the right-hand column. In this befiment, the chamber was inoculated 6.5 hours after the gradient was initiated. As in n l, chemotactic migration toward the aspartate source and formation of a chemotactic l & five were evident in both the experiment and the model predictions. However, the T es of the patterns were significantly different in the two runs. In particular, the \‘\] “ *wth pattern for Run 2 was more elongated and exhibited less flattening as it iroached the source. Figure 'l 38 . ma‘wmmmmwwmmm‘mmmwx Figure 13. Simulation compared to experiment for 1.0 mM aspartate source 39 Figure 14. Simulation compared to experiment for 0.1 mM aspartate source. 42.7 Effects of The simu equally spaced n: of the interval si and Run 1 was s two central axes on the predicted did show a mark tally .‘l [(D ‘ (u/rc‘) ( ‘0" do- I:igure 15 4'2'8 Content The ce‘ COUCemraUOH ] All fOur gram 40 Effects of grid-spacing on simulations The simulations shown in Figure 13 and Figure 14 used a grid consisting of 81 ly spaced nodes in each direction, for a total of 812 (6561) points. To test the effects : interval size on the predicted profiles, the grid was reduced to 512 (2601) points, Lun 1 was simulated again. Figure 15, which gives the cell concentrations along the entral axes of the DGC, shows that the increase from 512 to 812 had almost no effect c predicted profiles. Other simulations (not shown) at grid-spacings of 252 and 152 row a marked difference from the results at 512 intervals. 2 _ O 51 intervals 7 T . 81 intervals 0 51 intervals 6 _ 2': 6 ._ -—81 inten/als U 39 5 ~- 2 T 3 4 -- N 5‘ ._ 8 “r .8 3 Q :2 23 a 3 4 U 1 __ 0 I r . -::==:::=::n O :3 0 1 2 3 4 5 0 1 2 3 4 5 x-dimension (cm) y-dimenslon (cm) Figure 15. Cell concentration profiles across the two central axes of the DGC. Concentration profiles of chemoattractants and nutrient The cell patterns developed in response to the underlying, time-dependent :ntration profiles of the chemoattractants and nutrient. Examples of these latter es are shown, along with the corresponding cell profile, in Figure 16 and Figure 17. our graphs have the concentration, g/cm3, on the vertical axis, and the spatial 4.2.7 Effects of The simu equally spaced n of the interval sf and Run 1 was 5 mo central axes on the predicted did Show a marl ‘ (ti/1T) os~ H density x I" 1 0.4 — v V Figure l 4'2'8 COncen The c concentration profiles are sf All fOUr grill 40 l.2.7 Effects of grid-spacing on simulations The simulations shown in Figure 13 and Figure 14 used a grid consisting of 81 :qually spaced nodes in each direction, for a total of 812 (6561) points. To test the effects >f the interval size on the predicted profiles, the grid was reduced to 512 (2601) points, md Run 1 was simulated again. Figure 15, which gives the cell concentrations along the .wo central axes of the DGC, shows that the increase from 512 to 812 had almost no effect )n the predicted profiles. Other simulations (not shown) at grid-spacings of 252 and 152 lid show a marked difference from the results at 512 intervals. 2 " o 51 intervals 7 " . O 51 intervals A _81 Intervals "‘ 6 __ ——-81 intervals 8 1.6 “ 8 5° 3 5 —~ 2 1.2 -- 2 4 _- ’i N 3* ,3‘ __ 'g 0.8 g 3 “g 1: 2 a- E 04 E! U ‘ u 1-- O , l o ,____"_"__,__ --- 0 1 2 3 4 5 0 1 2 3 4 5 x-dimenslon (cm) y-dimension (cm) Figure 15. Cell concentration profiles across the two central axes of the DGC. 4.2.8 Concentration profiles of chemoattractants and nutrient The cell patterns developed in response to the underlying, time-dependent :oncentration profiles of the chemoattractants and nutrient. Examples of these latter profiles are shown, along with the corresponding cell profile, in Figure 16 and Figure 17. All four graphs have the concentration, g/cm3, on the vertical axis, and the spatial position. in cm. graphing prograi 41 sition, in cm, on the horizontal axes. Note that the jagged edges are artifacts of the aphing program. 'F‘ 42 Figure 16. Top: Typical cell profile. Bottom: nutrient profile. 43 “was; ‘ “was“ ‘smm‘“““ “‘ ‘\“\‘ \ o \‘ .‘¢ ‘ . “\‘3 ‘ “t\“““ ‘ “\““‘ "“‘.“ .o a,“ ,a Figure 17 To ' p. T ‘ ypical S profile. Bottom : Q profi 1e. Experimental VI Figure 18 (’Wid oxygen profiles away from the 2 Oxygen Concentration (p M) T0 val measured in a the Same as . compared to ; [CSUhS 0f the 44 imental verification of these profiles was more challenging. However, as shown in 3 18 (W idman et al., 1997), microelectrodes have been used to confirm that sharp :n profiles exist across both the cell from moving toward the aspartate source and from the aspartate source. 250 lnsrde cell population pl A a .. -o I .d A I a . r V I C I '2 ID ' e 150 "" I l 5 55' . 8 . o A o 100 -~ I / c l a, I a) ’ 5 Q 9 / o I A 50 « (I ’ ‘ I 4'- -A- 4 ’ Outside cell population 0 r r i i . -6 -4 -2 o 2 4 6 Distance from Wave Front (mm) Figure 18. Oxygen gradients in the DGC arena To validate the model's predictions of chemical gradients, glucose gradients were .1de in a miniature DGC using microbiosensors. The miniature DGC is fashioned me as a regular DGC, but has an arena length and width of only 3.0 cm, as ared to 5.0 cm for the regular DGC. Figure 19 (Widman et al., 1997) shows the s of these measurements for a single time point, along four different lines in the miniature DGC. openings not C( model accounts 45 ;e DGC. The variations in the profiles for the different lines arise from reservoir :3 not completely spanning the width of the DGC, as shown in Figure 12. The ‘zcounts for this geometry. 257. n U a I I n.- I a I :2...- . 0.500 OIOU-a-mv ‘ A i u... ”w I)» 1 v 1 :21. -UCOU avg-003.0 46 Line 1: 1.5 mm from top in 5 10 15 20 25 I Position [mm] Line 2: 7.5 mm from top 2500 2000- 1500-- 1000-- 500-- Glucose Conc. firM] Position [mm] Line 3: 12.5 mm from top L r r a 4' o 5 1o 15 20 25 Position [mm] Line 4: 17.5 mm from top 2500 2000‘- 1500-~ 1000- 500~~ Glucose Conc. [ng] Position [mm] Line 2 Line 1 Glucose Source -—-~-~-—-~44-----~ m \\ Line 3 Line 4 / 4' Figure 19. Glucose gradients at several positions in the DGC 4.2.9 Assumpt The mor and boundary c between the er assumptions USI In the l depend upon b neglect these c random motilit to chemotaxis. relatively inser Anothe bacteria are n Equation ( 13)) assay (Ford e1 C0ncentrations Chemoattractar prevents the 11 Cells. In Our r the Simlllatior speed of 7.92 gradient 51ml: 47 ptions and simplifications in the mathematical model '3 model was able to reproduce several experimental trends for different initial ' ary conditions using the same set of parameters. The reasonable agreement 1e experimental and modeling results provides support for the simplifying : 3 used. These assumptions are discussed below in more detail. he RTBL model of chemotaxis, the random motility coefficient is shown to ;on both temporal and spatial chemoattractant gradients. We have chosen to ese dependencies and assume constant 11. The magnitude of the flux due to )tility is predicted to be several orders of magnitude smaller than the fluxes due lXiS. In results not shown, the predicted growth patterns have been found to be nsensitive to it. Therefore an exact calculation of u is not deemed necessary. other assumption is that the gradients of chemoattractant encountered by the re relatively shallow (i.e. Equation (14) is a reasonable approximation to 13)). This assumption has been found to be adequate for modeling the capillary ‘d et al., 1990; Rivero and Lauffenburger, 1986). However, the locally high ions of cells in the chemotactic waves are associated with significant ictant gradients. The use of the hyperbolic tangent term in Equation (13) 1e predicted chemotactic velocity from exceeding the swimming speed of the hr modeling simulations, the maximum predicted velocity at any time during .tion is only 0.154 cm/h (0.43 uni/s), much less than a typical cell swimming :92 cm/h (22 urn/s) (Frymier, et al., 1994). To further validate the shallow- ,implification, velocity profiles from the hyperbolic tangent model and the shallow gradien gradients of S 1 (Ford and laufl swimming spec the number of i used in the sim sensitivity of 6. in Figure 20 to introduced by encountered in 48 ow gradient model were calculated by Equations (13) and (14), respectively, for the ients of S predicted by the model. Typical values for the individual cell parameters :1 and Lauffenburger, 1990) for E. coli were chosen for use in Equation (13). The cell .nning speed, 1), was taken to be 30 urn/s, and the tumbling frequency multiplied by pumber of receptors, ONT, was 75 s. The value for KD, 2.0x10'6, was the same as that :in the simulations. Equation (15) was used to calculate a value for the chemotactic tivity of 6.74x10'4 cm2/s from the single cell parameters. The comparisons are shown gure 20 for the first and last time points in Runs 1 and 2. There is virtually no error duced by substituting Equation (14) for Equation (13) for any of the gradients untered in either run. (1‘ (3}). . (.‘i ((7: U (o :2 Eu: -G\ Eer-U «undo « >«.00_0> 0. 0 0. w ~I>t Figure 49 7.0E-05 -- Run 1 —Tanh model 6.0E-05 -l 0 Shallow gradient model 5.0E-05 - l 4.0E-05 - l 3.0E-05 - 2.0E-05 - Chemotactic velocity to S [cm/s] 1.0E-05 - ctctcttccg, .‘ 0.0E+00 Jagauuuuu 'c. i 1 -1 .0E-05 y-Positlon [cm] 7.0E-05 -— Run 2 6.0E—05 -— 5.0E-05 _. 0 Shallow gradient model 4.0E-05 - 3.0E-05 a 2.0E-05 - Chemotactic velocity to S [cm/s] 1.0E-05 — 0.0E+00 ‘ " . l -1 .OE-05 y-Posltlon [cm] Figure 20. Comparison of model with and without hyperbolic tangent term. A third constant throu uniformly dist gradients to f significant ve rnicroelectrode partially respo: 4.3 Analysi 4.3.1 Flux ca The m from chemota 0f fluxes due t Where J W is th [0 S'gradients magnllllde am 20 hr and are auow COITesI ConesPonds n 50 A third assumption is that the concentration of all components in the model is instant throughout the depth (the z-direction) of the gel. Oxygen, which is initially 1iformly dispersed in the gel, is quickly consumed by the cells, causing vertical "adients to form. In experiments (Emerson, unpublished data), we have measured gnificant vertical oxygen gradients in the zone of cell growth with oxygen icroelectrodes. The variation in oxygen concentration with depth is thought to be utially responsible for the deviations between the experimental and predicted profiles. .3 Analysis of modeling results 3.1 Flux calculations The model was used to calculate the direction and magnitude of cell fluxes arising om chemotaxis and random motility in the DGC. Equation (11) was rewritten in terms i'fluxes due to single driving forces as Ju : Jim + JuS + JuQ (26) here J W is the flux due to random motility, J uS is the flux due to chemotaxis in response S-gradients, and J llQ is the flux due to chemotaxis in response to Q-gradients. The agnitude and direction of each flux term were calculated individually for Run 1 at time ) hr and are presented at each node point in Figure 21—Figure 24. The length of each row corresponds to the magnitude of the flux, and the direction of the arrow >rresponds to the direction of the flux. The tail of the arrow 1s located at the pomt where e flux is calculated. (Note that in Figure 21-Figure 24 only the top 20% of the fluxes are shown to reduc the lighter (wh illustrated. Th chemoattractan arrows onto th profile. and th Figure 33 sho' obtained by ad type of cell tlu 51 vn to reduce the number of arrows and to improve clarity.) In each of these figures, ighter (whiter) shading corresponds to higher concentration of the component being trated. The top graph in Figure 21 overlays the J uS flux arrows onto the moattractant S gradient profile, and the lower graph in Figure 21 overlays the J us flux ws onto the cell profile. The top graph in Figure 22 shows J uQ overlayed onto the Q ile, and the bottom graph in Figure 22 shows J uQ overlayed onto the cell profile. re 23 shows Juu, the flux due to random motility. Figure 24 shows the total flux, ined by adding the three flux components together. The maximum magnitude of each of cell flux is given in Table 3. Table 3. Maximum predicted cell fluxes for Run 1 at 20 hours. Flux component Maximum magnitude of flux (x 105 gcells/cmz‘h) Jus 8.05 J uQ 3.83 J m, 0.0000131 J u 1 1.4 52 5 2 FE: conical» x—position [mm] 5O 5 2 7:5 coEonl> x-position [mm] Figure 21. Fluxes due to response to S chemoattractant. y—position [mm] 25 x—position [mm] 50 y-position [mm] N 01 25 x—position [mm] Figure 22. Fluxes due to response to Q chemoattractant. y—position [mm] N 01 25 x—positlon [mm] Figure 23. Fluxes due to random motility. y—position [mm] l\) U! 25 x—position [mm] Figure 24. Total mass fluxes. 432 Receptl The (16 be calculated 1 05. as the aver Sdivided by t TECCplOTS WCTI The value of response is c zero. the chei factor has be 13 and Figu Calculation ( of the maxirr 55 3.2 Receptor saturation The degree to which receptor saturation suppresses the chemotactic response can : calculated from the simulation results. We define a global chemotactic response factor, ;, as the average cell flux (along a line of constant x) in response to the chemoattractant divided by the average cell flux (along the same constant x line) that would result if the :ceptors were completely free of any S. This definition is expressed mathematically as 2R 2R 2 ‘2 (27) KB, 35 R s as {[9503 (K03 +S)2 3y “idy lid); iii“- KDS] 3y uddy 5 = 2R 2R : 2R 1 8S 3S ———u d d —u d {[9605 KDSBy :iy ii y {[3}; iy he value of ¢s will vary between 0 and 1. If tbs is close to one, then the chemotactic spouse is close to the maximum it can attain for the given gradient. As (1)3 approaches arc, the chemotactic response to the gradient of S diminishes. The chemotactic response lctor has been calculated for Run 1 and Run 2 at each of the time points shown in Figure 3 and Figure 14. Figure 25 shows (pg plotted against the x-position for both runs. alculation of (pg was limited to areas where the cell concentration was greater than 0.1% Fthe maximum cell concentration occurring at that time point. 0.2 4’s Figure 2 56 0.20 ' ' ' 11 hr —0— 17 hr 5 —20 hr 0 N I ‘ l \‘ ,I \ —O—23 hr V B B (D e 0.10 - ‘_ ‘_ 0.00 ‘r .1- 0.90 I - - - 9hr —0—15hr —18 hr —O--24 hr I 0.80 -— ' 0.70 ~— 0.60 -- 0.50 -_ s 0.40 -+- 0.30 -- 0.20 *- O.10 ~F 0.00 —h— ‘1— .- x-position [cm] Figure 25 Top: $3 for Run 1. Bottom: (pg for Run 2. Note different vertical scale. 4.3.3 Pattern Twoh first hypothCS reducing the c takes longer f near the sourc The cl first hypothes 03 was deterr point A at tin of the patterr values vary 1 lesser degree the maximun Run l. The 2 While the 24 occur, the y: in the center. Suggest that ; 57 4.3.3 Pattern dynamics Two hypotheses could explain the flattening of the pattern seen in Figure 13. The first hypothesis is that the cell receptors become saturated for the chemoattractant S, reducing the cell’s chemotactic response to gradients in S. The second hypothesis is that it takes longer for the cells to consume the higher concentrations of chemoattractant found near the source. The chemotactic response factor values (Figure 25) provide strong support for the first hypothesis. In Run 1, (p3 values lie between 0.02 and 0.18 at all time points for which (1)3 was determined. The values are the lowest in the middle of the pattern (indicated by point A at time 17 hr), where the flattening is the most evident, and highest near the edges of the pattern (indicated by the two points labeled B at time 17 hr). In Run 2, the tbs values vary between 0.25 and 0.9, indicating that chemotaxis is suppressed to a much lesser degree by receptor saturation than in Run 1. This result would be expected since the maximum chemoattractant concentration in Run 2 is approximately one tenth that in Run 1. The 23 hr (pg curve for for Run 1 is flat from approximately 1.7 cm to 3.2 cm, while the 24 hr (1)3 curve for Run 2 has a much rounder appearance. For broadening to occur, the y-component of the wave’s velocity at the sides of the pattern must exceed that in the center. The shapes of the $3 curves are consistent with such a velocity gradient, and Suggest that receptor saturation is at least partially responsible for the pattern broadening. 4.3.4 Chemo Mover source reservt was not reprt that this beha second chemt second chem coli (Adler. 1 significant or away from th moving away gradient. 44 Variat The exllellments Compounds, however’ an ComPlelt m0 First COmpetition modificatim 58 4.3.4 Chemotactic wave in response to oxygen Movement of a portion of the chemotactic wave away from the chemoattractant source reservoir (i.e. toward the north) was an unexpected experimental result. This effect was not reproduced in the simulations with a single chemoattractant. We hypothesized that this behavior could arise in response to gradients formed by the consumption of a second chemoattractant that is uniformly distributed throughout the chamber initially. The second chemoattractant has been identified as oxygen, a known chemoattractant for E. coli (Adler, 1972). The microelectrode measurements shown in Figure 18 indicate that significant oxygen gradients coincide with the cell wave fronts moving both toward and away from the applied chemoattractant gradient. Model predictions indicate that the wave moving away from the aspartate source tracks an oxygen gradient rather than an aspartate gradient. 4.4 Variations on original model The original model, presented in Section 4.2, gave good agreement with experiments for one bacterial population responding to gradients of two chemical compounds, and growing on a third. The resulting simulations were not perfect matches, however, and experimental data, such as the oxygen gradients, suggested that a more complex model was needed to match the experimental results more closely. First, a second bacterial balance was added to the model. This balance allowed competition or other multi-population phenomena to be studied. More details on this modification to the model are given in Sections 5 and 6. In addi second-genera and the nutrie speaking, this Q were oxyge oxygen. If ( improvement and 01115. 19E takes into 3C1 period of it chemoattract; growth. The motility. for. All it cells to grol comPound u coli, or acetz The : du. we Where Vei l popuiaUOn‘ 59 In addition to the second cell balance, several improvements were built into this second-generation model. These included making consumption of the chemoattractant S and the nutrient H dependent on the concentration of chemoattractant, Q. Biologically speaking, this can be interpreted as modeling Q as the electron acceptor in the system. If Q were oxygen, then the cell's metabolism is aerobic, with electrons being donated to oxygen. If Q were nitrate, anaerobic denitrification could be modeled. Another improvement was the addition of death terms and endogenous metabolism terms (Bailey and Ollis, 1986). The death terms model the death of cells, while endogenous metabolism 1 :takes into account the loss of cell mass due to internal consumption, such as during a period of low nutrient availability. Cell maintenance terms were added to the chemoattractant and nutrient balances to allow for consumption for uses other than growth. The cell maintenance terms could include nutrient consumed for energy for motility, for example. An important modification to the second-generation model was the ability of the cells to grow on the chemoattractant S. This gives the model the ability to model a compound which is both a chemoattractant and a growth nutrient, such as glucose for E. coli, or acetate for Pseudomonas KC (see Section 7). The second-generation model for cell population i (i=a or b) is given by 8L1, 797:“:‘V2 up _V.(ch)is _V.(V It)“, + 1121‘ i . . 2 Q vth + vlSS __ veg u,- __ d Ll ( 8) C,,, + H C“ + S where vet is the endogenous metabolism term and dth; is the death term for the ith population. The balance for the consumable compound j (i=8 or H) is now given by where vpm is term is given 1111616 ij IS 1 term. for the 1 (9Q _ E _ u‘here Yiq is endogenous 1 Two ShOWn in p metabolism 1 to Obtain ac C 60 if 2 1' Q (29) — V 31—11) j 2["”'C+JC +Qu'] where Vijm is the total consumption coefficient for the ith population growing on j. This term is given by v.. ( 30) ijm — Y ij '1 where Yij is the yield coefficient for cells growing on j and cmij is the cell maintenance term, for the consumption of j. The balance for the compound Q is given by mg a Q _ rhm H vism S vet Q ( 31) BI—qu—VQ 21“th C,,,+H+Y C+S+YiC,q+Q“‘] where Yiq is the yield coefficient for the consumption of Q by the ith population used for endogenous metabolism. Two time-points from a sample simulation using the second-generation model are shown in Figure 26. Values for unknown parameters, such as the endogeneous metabolism terms and the consumption terms were estimated. More work will be needed to obtain accurate values for these parameters, but the model has been shown to produce results that resemble the experimental data. 4.5 Laser The 1 this research Student, Seb; Will be presl 1997. Moti COtifficiem : moving thro Therefom, 3 allows the r: 61 Figure 26. Sample simulations from second-generation model. 4.5 Laser diffraction capillary assay The laser diffraction capillary assay was conceived and developed as a portion of this research project, but the majority of the work was done by an Aachen exchange student, Sebastian Schmidt. For this reason, only a brief summary of the LDCA method will be presented within the body of this work. For further details, see Schmidt et al., 1997. Motility parameters necessary for the modeling work, such as the random motility coefficient and the chemotactic sensitivity coefficient, were not known for bacteria moving through a semi—solid medium, such as the dilute agar gel in the arena of the DGC. Therefore, a method to independently measure these parameters was needed. The LDCA allows the random motility coefficient to be measured at various agar gel concentrations, mdmmmm coefficients It Thel mounted on 1 was achieved uupenuonsi mwmmnh alle/Ne lase experiment. 62 and preliminary experiments indicate that it will also allow chemotactic sensitivity coefficients to be measured as a function of agar concentration. The LDCA consists of a 1.6 mm inside diameter glass capillary, which can be mounted on a microscope stage. A step change in cell concentration inside the capillary was achieved by successively inserting the tip of the capillary into each of the two agar suspensions in the culture tubes (Figure 27). The capillaries were filled over a total length of 30 m, then mounted on the microscope stage. The end of the capillary opposite from a He/Ne laser was sealed with silicone grease to avoid unwanted convection during the experiment. The laser w; Was mounte 15 minutes, mElintained ; The Latlffenburg dimenSlOna] 63 x=-L ' x=0 =L = Agar suspension with no cells - = Agar suspension with cells Figure 27. Filling the capillary in the LDCA The laser was shone axially through the capillary. A CCD camera, attached to a computer, was mounted on the optical lens of the microscope. Images were captured at 1.5 minutes, 15 minutes, and then at 30 minute intervals. A typical run lasted 3 hours. The LDCA was maintained at room temperature for all runs. The mathematical model for the LDCA is similar to that used by Ford and Lauffenburger (1992) for the SFDC, with the chemotaxis terms omitted. The one- dimensional cell balance is given by The analytica using the b0t The concem This SiOpe c The model i m01ility We a function ( Points, The confidence 64 a_B_ 323 (32) at flaxz The analytical solution to this model is B 1 .. (33> ——2-erfc Tilt using the boundary conditions B(x,0)=1for—L<10'6 cm2/s). The important finding of the LDCA work was that the random motility coefficient displayed a linear relationship with agar concentration, at least in the region of agar concentration tested. 25 11 x107 (cmzls) 0 h r 0.1 0.2 0.3 % Agar Concentration Figure 28. The random motility coefficient of Pseudomas KC. 5. MICRC 5.1 Compl Cherr to gradients 1 119761 Stated Since most e sole criterion Thee utilizing mat microbial pt confined do) motility Cou 01- (1982') w Of chemotax growth rates bacterial p0) (Lauffenbm populations Sufficiently had a highe geometry ‘ 5. MICROBIAL COMPETITION 5.1 Competition background Chemotactic microorganisms able to position themselves optimally with respect to gradients may have a competitive advantage over other organisms. Chet and Mitchell (1976) stated that chemotaxis "presumably gives. . .microorganisms a selective advantage. Since most ecosystems are not fully mixed, enzyme kinetics cannot be postulated as the sole criterion for the competitive advantage of one microorganism over another." The effects of chemotaxis on competition have been studied on a theoretical basis, utilizing mathematical models to predict possible outcomes of competition between two microbial populations. Lauffenburger et al. (1981) showed that in a one-dimensional, confined domain, with a substrate diffusing in from one side, increasing levels of random motility could be detrimental to the survival of a bacterial pOpulation. Lauffenburger et al. (1982) went on to show that if the substrate was also a chemoattractant, a certain level of chemotaxis could impart an advantage to a population, canceling the effects of lower growth rates or higher random motilities. Both of these studies involved only a single bacterial population. Lauffenburger and Calcagno (1983), building on the previous work (Lauffenburger et al., 1981; Lauffenburger et al., 1982) modeled two randomly motile populations simultaneously. They found that if a slower growing population also had sufficiently lower random motility, it could out-compete a faster growing population that had a higher random motility. Kelly et al. (1988), using a model systemiwith the same geometry as the Lauffenburger papers, modeled two chemotactic pOpulations 66 simultaneous] above which immotile pop The n steady-state : first populati would survil authors poin order of a ye occur over 1 motility ant environmen Rela from an ex] 51min of Pr CUlture, bur (1985) stud Strain of 11. Chemotacti Similar gro Otl Organisms 67 simultaneously. Their model predicted that there was a minimum level of chemotaxis above which a population would have a competitive advantage over a non-chemotactic, immotile population. The model of Kelly et al., which included a death term and therefore allowed for steady-state solutions, predicted that at steady-state, three conditions could arise: (1) the first pOpulation would survive and the second would disappear, (2) the second population would survive, while the first disappeared, or (3) both populations would coexist. The authors pointed out that the time to attain steady-state predicted by their model is on the order of a year, and in real systems, it is likely that variations in chemical gradients would occur over that time period. They concluded that the ability to "consider the effects of motility and chemotaxis on population growth and competition in rapidly changing environments is clear". Relatively few studies have addressed the effects of chemotaxis on competition from an experimental vieWpoint. Pilgram and Williams (1976) found that a chemotactic strain of Proteus mirabilis outgrew a non-chemotactic but still motile strain in stationary culture, but that the two strains grew equally well in mixed’culture. Kennedy and Lawless (1985) studied a motile, chemotactic strain of Pseudomonas fluorescens and a immotile strain of the same species, and found that in unmixed aerobic and anaerobic soils, the chemotactic strain survived significantly better than the non-chemotactic strain. Both had similar growth characteristics in mixed culture. Other work has focused on the interaction of two or more populations of Organisms growing in the presence of spatially varying environments. Caldwell and Hirsch (1972 gradients int layer. and grl study did no could have l the local con Tilm. structured he dispersal rat movement c studied how competition In t1 Spatially 51p. Called the d the format Microorgan movement mathematic balance, an Widman et (1) the $111): 68 Hirsch (1972) studied the growth of microorganisms in two-dimensional concentration gradients in their steady-state diffusion plate. The organisms were immobilized in an agar layer, and growth was quantified as a function of position in the gradients. Although this study did not address competition directly, the researchers did find that the organisms could have vastly different growth characteristics, depending upon their adaptability to the local concentration conditions. Tilman (1994), taking a broad ecological perspective of competition in spatially structured habitats, noted that coexistence of species could occur if a species with a high dispersal rate could move into a zone not inhabited by a superior competitor with less movement capability. Holmes et al. (1994), again from an ecologically based viewpoint, studied how partial differential equations could be used to study, among other things, competition, dispersal, and dispersal-mediated coexistence in spatially structured systems. In this study, a system of partial differential equations was used to describe a spatially structured system in which two bacterial populations were growing. The system, called the diffusion gradient chamber (DGC), has been established as a tool that allows the formation of two-dimensional chemical gradients (Emerson et al., 1994). Microorganisms inoculated into the chamber can be observed, and their growth and movement properties recorded. Widman et al., (1997) (see Section 4.2), developed a mathematical model of the DGC that included two chemoattractant balances, a nutrient balance, and a cell balance. In addition to modeling a different physical geometry, the Widman et al. model offered three features not found in the Kelly et al., (1988) model: (1) the substrate does not necessarily have to be a chemoattractant, and vice versa; (2) the model equati transient m0 analyzed wit extended to competitions 5.2 Comp The system. intr balance. Th- Where, for i is the nutri Coefficient the receptc KDiQ is the gTOWth rat. Th. the second 69 model equations are solved in two dimensions, as opposed to one dimension; and (3) the transient model can be solved over relatively short time periods, and the data can be analyzed without requiring a steady-state solution. In this study, the model of the DGC is extended to include a second cell balance, which allows predictions of the outcomes of competitions between two microbial populations. 5.2 Competition mathematical model The system of coupled partial differential equations that describes the DGC system, introduced in Widman et al. (1997), was extended to include a second cell balance. The cell balances for Populations A and B are given by 38 W ( > au. K . KDiQ —'=fl.-V2u,-_Z.~ V __‘L’ uiVS —Zi V _— uin + ”i at OS [[(Kots +S)2 OQ (KDrQ +Q)2 C1” + H where, for i=A or B, u, is the cell concentration, S and Q are the two chemoattractants, H is the nutrient, u, is the random motility coefficient, X015 is the chemotactic sensitivity coefficient to the attractant S for the ith population, KDiS is the dissociation constant for the receptor-S complex, inQ is the chemotactic sensitivity coefficient to the attractant Q, KDiQ is the dissociation constant for the receptor-S complex, V,“ is the maximum specific growth rate, and CiH is the half-saturation constant for growth on H. The chemoattractant and nutrient balances were changed to reflect the addition of the second cell population. The chemoattractant balances are given by where, for j: chemoattract chemoattract chemoattract As (1 more of the chemoattrac the reservoir Which is a cl The Where DH i glOWIh of tl The batted in balance are nutrient are “11 popula where, for j=S or Q, Dj is the chemoattractant diffusion coefficient, Vij is the specific chemoattractant consumption coefficient for the ith population consuming the jth chemoattractant, and Cij is the saturation constants for ith population consuming the jth chemoattractant. As described in Section 4.2.3, the chemoattractant S is introduced from one or more of the DGC reservoirs, and forms diffusion gradient across the gel in the arena. The chemoattractant Q is initially present at a constant concentration throughout the gel and the reservoirs. In the experiments used to validate the model, Q was most likely oxygen, which is a chemoattractant for many bacteria. The nutrient balance becomes flzDHVZH—X B VrHH _L_‘i_ (40) at i=A C... + H Y... where DH is the nutrient diffusion coefficient, and YiH is the yield coefficient for the growth of the ith population on H. The assumptions made in the original, single population model are given and justified in Section 4.2.9. The new assumptions introduced by including the second cell balance are that the effects of the two populations on the chemoattractants and the nutrient are additive (Lauffenburger and Calcagno, 1983; Kelly et al., 1988); that neither cell population preferentially consumes one chemoattractant or nutrient over the other; and that the: predation or 5.3 Mode? .Vlanj parameters 1 correlations. model was i In th 71 and that there are no direct interaction effects between the two populations, such as predation or parasitism. 5.3 Modeling parameters Many of the parameters used in the model are given in Section 4.2.5. The parameters were measured independently, obtained from the literature, or calculated from correlations. All the parameters had reasonable values for a population of E. coli, and the model was verified by comparing its simulations to experiments with E. coli. In this study, the base set of parameters for each population is given in Table 4. Frt their effec Where the parameter 72 Table 4. Base case parameters for competition simulations. Parameter Value (i = a or b) u, 0.0010 cm2 h'1 XOAQ 0.085 cm2 h“ XOBQ 0 XOAs 0.020 cm?1 Xoas 0 Cir; 4.08x10’6 gH cm'3 CiQ 6.70x10'8 gQ cm’r Cis 5.50x10'5 gs cm‘f DH 0.01 c? h'1 DQ 0.033 cmTh’l D5 0.033 cmThT KDtQ 3.301.105 g., cm'3 Kors 2.00x10'vgs cm'3 vAH 0.35 h'1 VBH 0.5 h'I VAQ 0'02 8Q gu'l h-r vBQ 0 MS 0.60 gs a.‘1 h'1 vBs 0 YiH 0.50 gu/gH From this base set, individual parameters for Population A were varied to test their effects on the outcomes of competition simulations. In many mathematical models, Where the effects of parameters are unknown, a dimensionless model enables many parameter effects to be studied at one time. In the case of this particular model, however, a dimensionless model did not decrease the number of parameters, and was therefore not useful for pa the attractant the specific chosen as 11 inoculation j simulations. while that it rate for P0p (Hansen ant including 0 that it is onl In 2 inoculation The base 0 Chamber, a 54 Anal 5.4.1 The Ty populallor Populatior [631111 is 73 useful for parametric studies (see Appendix G). The chemotactic sensitivity coefficient to the attractant S (xOAs), the dissociation constant for the receptor-S complex (KDAS), and the specific chemoattractant consumption coefficient for the attractant S (vAs) were chosen as the variables to be manipulated. The initial amount of H present and the inoculation pattern of the two populations were also varied in some simulations. For all simulations, the maximum specific growth rate of P0pulation B (ttB) was set to 0.5 h'l, while that for Population A was left at the base case value of 0.35 h]. The larger growth rate for Population B would give it a competitive advantage in a well-mixed environment (Hansen and Hubbell, 1980). Population B does not interact with either chemoattractant, including chemotaxis or consumption. In other words, the base case for Population B is that it is only able to move by random motility, and it consumes only H. In all simulations, the gradient was allowed to initialize for 6 hours before inoculation. After inoculation, the duration of the simulated experiment was 30 hours. The base case initial condition for the cells is that they are inoculated in the center of the chamber, at an equal concentration, as described in Widman et al. (1997). 5.4 Analysis of competition simulations 5.4.1 The dynamic competition factor Typically, the result of a competition experiment between two microbial pOpulations in well-mixed culture is quantified by the death or disappearance of one population and the continued survival of the other. This type of steady-state competition result is not readily amenable to the transient DGC model. Therefore, a dynamic competition masses of th where V is ‘ the model. dimensional depth of the W1. Value of W a Value of 74 competition factor (tlt), was defined to give a time-variant measure of the ratio of the masses of the populations; ‘1’ is defined by ”I“ (x, y,z,l‘)dV (41) 1”“) = yyyu,(x,y,z,.)av V where V is the total volume of the gel, and ua(x,y,z,t) and ub(x,y,z,t) are calculated from the model. The coordinate axes are defined in Figure 29. It is assumed that the two- dimensional concentration profiles predicted by the model are constant throughout the depth of the gel (the z-direction). Source side X Sink side Figure 29. Coordinate axes used for competition analysis. When \p = 1, both populations have the same total mass present in the chamber. A value of tyr < 1, indicates that Population B has a higher total mass than Population A, and a value of \p > 1 indicates that Population B has a lower total mass than Population A. Therefore, 1 instantaneou time. and vit Anot derived fror masses. and lftp‘>0,tl B.lftp<1 although its 5.4.2 The Ant P0pulation average cc X-directior COllcentrat 75 Therefore, for the purposes of this study, Population A will be said to have an instantaneous competitive advantage over Population B at a given time if \p > 1 at that time, and vice versa. Another measure of the competitive relationship of the two populations can be derived from Ill. The local competition rate, ty‘, is the rate of change of the ratio of the masses, and can be found by , d '1’ (42) If 111' > 0, then the mass of Population A is increasing relative to the mass of Population B. If t)! < 1, but 111' > 1, Population A may still have the overall competitive advantage, although its total mass has not yet become greater than the total mass of B. 5.4.2 The average concentration curve Another useful approach to studying the spatial arrangement of the two competing pOpulations was the development of the average concentration curve. This curve plots the average concentration of the population along each line perpendicular to the gradient (the x-direction) vs. the position parallel to the gradient (the y-direction). The average concentration is calculated (for Population A, for example) by x=L n n . (43) j uA(x, y)|yAyAzdx Ema, y)|yAyAzAx,. 2...,(1, y)) Y , i=1 : i=1 ' HA(Y) = x=0 AxAyAz z nAxAyAz n where L is t1 and Ay are . the gel. Figt point on the image. TOT any 1 ConCentra the Simul: 76 where L is the width of the chamber, n is the number of grid points in the x-direction, Ax and Ay are determined by the grid spacing and are always equal, and Az is the depth of the gel. Figure 30 shows a typical competition simulation, and one line used to calculate a point on the total mass curve. The chemoattractant S is diffusing in from the top of the image. 15 hours 20 hours 25 hours Line of constant y along which one point of average concentration curve calculated Figure 30. Typical simulation, illustrating average concentration curve calculation. It should be noted that the average concentration curve can be related to the dynamic competition factor as jit‘A(y,t)dy (44) WU): j'rZB(y,t)dy for any time point. By studying the dynamic competition factor and the average concentration curve, some measure of the spatially developing competitive outcomes of the simulations can be found. 5.4.3 Satur- The utilized to e: quantity use fora given 1 the equality Where Fm, of constant 5-5 Com Ma W33 also 1 ”We. and The Matl. appTOle; 55.1 an In in Table 4 77 5.4.3 Saturation of chemoreceptors The global chemotactic response factor, 03, defined in Section 4.3.2, can be utilized to explain some of the responses that occur in the following simulations. Another quantity used to analyze some of the simulation results was the maximum flux attainable for a given S gradient. This term is the integrand of the denominator of the second term of the equality in Equation ( 27), and can be defined as 10.45 é‘gd (45) y KDAS 8y x L W O F,,,,.(x) = where Fmax(x) is the average maximum flux attainable for a given S gradient, along a line of constant x. 5.5 Competition simulations Matlab v. 4.2c1 was used to visualize the output of the FORTRAN model. Matlab was also used to solve for the dynamic competition factor, the average concentration curve, and the global chemotactic response factor (see Appendix F for the Matlab code). The Matlab function "trapz.m", which utilizes the trapezoidal method, was used to approximate the solution to the integrals. 5.5.1 Effect of modeling parameters on competition results In the following sets of simulations, the parameters for Population B are as shown in Table 4, The indicated parameter for Population A is varied in each set of simulations. Range findi showed a gt range to prt concentratit hours after prior to cel number is 1 al.. 1997). simulation 5.5.2 Var: Tilt 0.01-0.071 points in Correspont gray. Pop Chemoattr the moon} dynamicc trend. 78 Range finding runs (results not shown) were performed to find a parameter range that showed a good variation in results. The parameter of interest was then varied within this range to produce the desired simulation sets. Both populations were inoculated in equal concentrations in the center of the chamber. The simulated experimental time was 30 hours after inoculation of the cells. The chemoattractant gradient was initiated for 6 hours prior to cell inoculation. The grid-spacing used for all simulations was 412 points. This number is less than the 512 points used in the study that validated the model (W idman er al., 1997), but was found to be sufficiently accurate and allowed for the fastest possible simulation times. 5.5.2 Variation of ans The chemotactic sensitivity factor for Population A, X0As. was varied between 0.01-0.07 cmzlh to assess its affect on the outcome of simulated competitions. Three time points in this simulation for X0AS=O~O7 are shown in Figure 30. The lighter areas correspond to higher cell densities. Both Populations A and B are shown in shades of gray. Population A is the one which moves out in a chemotactic wave toward the chemoattractant source, while P0pulation B remains in the center of the chamber, around the inoculation point. Figure 31A shows the dynamic competition factor. In many of the dynamic competition graphs shown, the initial 10 hour period shows a similar decreasing trend. Thi center of t] Figure 31E graph and the positic Chamber 1 reached a. limb Popu C0“Central so they W hour time have Th. the Chemt 79 A XOAs=0.07 —U— ___ . :XOAs=0.05 1-4 -- chs 0 O1 XOA5=0'03 ......... xOAs=0.05 XaAs=0-03 o XOAs=0-07 ”o ,_ X < ‘5 C .2 A "' (3‘ 5 08- g g o... 9‘ 5 a. 0.6 3 V 0.6 o 04 3. ° 3’ 0.4 0'2 ----I a 02 II... > ' 'I'llll... x0A8=O 01 < o r 4 o 0 10 20 30 40 Time (h) y-Position (cm) Figure 31 Variation of XOAS- This trend occurs because significant concentrations of S do not diffuse into the center of the chamber until approximately 10 hours after inoculation (results not shown). Figure 31B shows the average concentration curves for four values of X0As, at 25 h. In this graph and all subsequent average concentration curve figures, the attractant source is at the position y=5 cm. For higher levels of XOAs. Population A encountered the wall of the chamber before 30 hours. In the graphs of Figure 31, only times before the wall was reached are shown in the graphs. Figure 32A shows the average concentration curve for both Populations A and B, for the case where dis=0-03 cmzlh, at 25 hours. The average concentration curves for Population B in most of the other simulations were very similar, so they will not be shown in future graphs. Figure 32B shows the density maps of the 25 hour time-point of Populations A and B for XoAs=0~03 cmZ/h, illustrating the chemotactic wave. The wave is the high concentration band of cells of Population A moving toward the chemoattractant source. o o . l 2*? .4.__.444U , Q N 3 ' 3 cm ) Cell concentration x10 (9! _O .,}U, A 0 Note that the highes This does that it is a T1 the cell p. Figure 31 grOWth a. I)OIJUIatir (To/13201 highest 1. Popuiaill 8O 0.5 61‘ E A ° l B) 0.4 T PopulationA v o P) :l a Population B . , E 3‘ l t: 0.3 .9 s i-l l E s‘ E 0.2 ‘y _‘ 0 E o 8 o 0.1 ‘1 T.» o n o ..... y ..... 0 1 2 3 4 5 y-Position (cm) Figure 32. Other results of variation of XOAS Note that in this figure, and all figures of the density map type, the cell population with the highest concentration at any point in the chamber is the one that appears in the graph. This does not necessarily mean that the other population does not exist at this point, only that it is at a lower concentration. The chemotactic sensitivity coefficient is a measure of the magnitude with which the cell population migrates chemotactically in response to a chemoattractant gradient. In Figure 31A it is apparent that below a certain X0 value (in this case, X0AS=O~01 cmZ/h), the growth advantage of Population B is so large that, in the time-frame of the experiment, Population A is overgrown by B, as indicated by 111' < O at all times. At higher levels (X0A520.03 cmz/h), the sign of 111' becomes positive after about 10 hours, and at the highest levels of XOAs. 11! reaches a value greater than one, indicating that the total mass of Population A has become greater than the total mass of Population B. For this set of conditions, . total mass 0 The the chemot average cor side of the Population peak of P0 two populz of higher l with B. Th1 sensitivity This indie Chiming c increasing Cells mov cells may and there CaITying . adVantagr and the n' 8 1 conditions, chemotaxis can be said to impart a competitive advantage, as measured by the total mass of microorganisms present, to Population A. The average concentration curves of Figure 31B give another view of the effect of the chemotactic sensitivity. As X0As increases, two things occur. First, the peak of the average concentration curve shifts toward the attractant source, which is located at the side of the chamber corresponding to y=5 cm. The average concentration curve for Population B remains close to the center of the chamber, as illustrated in Figure 32A. The peak of Population A moving further from the center indicates greater separation of the two populations. By moving away from Population B, Population A is able to enter areas of higher nutrient (H) concentration, that can be consumed without having to compete with B. The second trend that can be observed in Figure 31B is that as the chemotactic sensitivity increases, the maximum height of the average concentration curve increases. This indicates that not only is the total mass of Population A increasing, but also that the carrying capacity of the chemotactic wave (Widman er al., 1997) of Population A is increasing. The chemotactic wave is better illustrated in Figure 32B. The bright band of cells moving toward the chemoattractant source is the chemotactic wave. This wave of cells may allow for faster consumption of the chemoattractant, creating sharper gradients, and therefore allowing for faster chemotaxis toward the source. By increasing the calTying capacity of the wave, the cells may be gaining an even larger competitive advantage by being able to consume larger concentrations of both the chemoattractant S and the nutrient H as the wave travels through the region. 5.5.3 Vari: The coefficient. dynamic cc cm) profile 33C shows 33D show 82 5.5.3 Variation of vas The second parameter varied was the specific chemoattractant consumption coefficient, vAs. It was varied between 0 and 1.0 gsguA'lh'l. Figure 33A shows the dynamic competition factor for five values of vAs. Figure 33B shows a centerline (x=2.5 cm) profile of the chemoattractant S gradient at 30 hours for each value of vAs. Figure 33C shows a density map profile for both populations with vAs=l .0 gsguA'lh'1 and Figure 33D shows the profiles for vAs=O gsguA‘lh'l, both at 30 hours. 0.1 ° 011 the d Flgm‘e ’ Other w. 83 1 1.2 0.9 — A B 0,8 - VAs=1-0 1 A 0.7 -- ”E 0.8 0.6 -- 0 .7: VAS=0.6 @ V 0.5 -- 0.6 9‘ ‘b 0.4 -' v- . VAs=0.3 X 0.4 lncreastng v.15 0'3 " VAs=0.06 ‘0 0.2 -- VAS=0 0.2 0.1 ‘r 0 ‘fi 1 l . 0 0 10 20 30 40 0 1 2 3 4 5 Time (h) y-Position (cm) Population A Population B Population A \ Population B Figure 33. Variation of vAs. Increasing the specific chemoattractant consumption coefficient, vAs, had an effect on the dynamic competition factor curves similar to that of increasing XOAs, as shown in Figure 33A. The reasons for the increased competitive advantage are not the same, however. Figure 33B shows the S concentration profiles for the position x=2.5 cm, or, in Other words, down the center of the DGC. As vAs is increased, a given mass of cells can consume : encountere response, vAg:0, me positive n were grea This can; of S, the When v, subseque: 33C, Whe and the ‘ consume (X0A3=0.( 5.5.4 v1 COmplex 343 she Show '1] 84 consume a larger quantity of S, and therefore increase the gradient of S that is encountered by the wave of A. This increase in slope enables a faster chemotactic response, and therefore a more rapid separation of P0pulations A and B. Even with vA5=O, meaning no consumption of the chemoattractant S, the value of w' does become positive near the end of the simulation (Figure 33A). If the dimensions of the chamber were greater, it is likely that w' would become positive, even with no consumption of S. This can again be explained by the carrying capacity of the wave. With no consumption of S, the gradient is relatively shallow, and the carrying capacity of the wave is low. When vAs is high, the gradients are sharper, which increases chemotaxis, and subsequently increases the carrying capacity of the wave. This is illustrated in Figure 33C, where VAS is high and the wave is very apparent, and in Figure 33D, where vAszO, and the wave is just beginning to appear at 30 hours. The ability of Population A to consume S gives it a greater competitive advantage, but at this level of chemotaxis (X0A3=O.02), the ability to consume S is not necessary for A to compete successfully. 5.5.4 Variation of KDAS The third parameter varied was the dissociation constant for the receptor-S complex, KDAS. The parameter range was 1><10'7 - 8><10'6 g/cm3. Figure 34A and Figure 34B show the w(t) values for five values of KDAS. Note that the curve for KDAS=1><10'6 is shown in both figures, for comparison. The arrows Show the direction of increasing KDAS. b-.. N.25 «:9 now x =5“. «asuxmv new x COSQbCOOCOO Onwfigw>< ‘P (t) Fun-x X 103 (90A ct“.2 hd) Average concentration x 10’ (gems) 1 A 0.8 u- 0.6 ‘ - F KDAs=1x10 6 0'4 -_ KDAhs=5X‘|0'7 0.2 _l Kms=mo" 0 i—‘ l r I 0 10 20 30 40 0.8 ‘- 0.6 < 0.4 0.2 ‘- O :3 05-. {Xx l. ) x-Position (cm) x-Position (cm) y-Position (cm) Figure 34. Variation of KDAS. Time=30 h for all figures. \< J ‘x «‘2! 0.1 )(xxxxxfl xxx X Xxxxxxx In ? competitiv 1x106, inc reasons fo average rr legend in the global maximum Consequc when S : points it decrease gradient. value fc Optimur. gradient the ave Value 0‘ Values interim 86 In Figure 34A, for KDAS < 1x106, increasing KDAS made Population A more competitive, as indicated by larger w values at the later time points. However, above 1x10'6, increasing KDAS caused \y to decrease at later times, as shown in Figure 34B. The reasons for these trends are shown in Figure 34C and Figure 34D. Figure 34C shows the average maximum attainable fluxes for the corresponding S-gradient, at 30 hours. The legend in Figure 34C applies to Figure 34D and Figure 34E, as well. Figure 34D shows the global chemotactic response factor for the five parameter values, also at 30 hours. The maximum mass flux is inversely proportional to KDAS, as shown in Equation (45). Consequently, at the lowest values of KDAS, Fmax reaches its highest levels. However, when S >> KDAS, (e.g. for the lowest values of KDAS), the value of (1)3 approaches 0 at all points in the chamber. As KDAS increases, the level of saturation of the receptors decreases, and therefore the population responds with a higher chemotactic flux to the S- gradient. The trade—off between the opposing trends of Fmax and (pg results in an optimum value for KDAS of about 1x10‘6 g/cm3, where flux is high and saturation is low. This optimum KDAs value is most likely dependent on the concentration and SIOpe of the gradient. Evidence for an optimum KDAS value is also given in Figure 34E, which shows the average concentration profiles for Population A for a low, intermediate, and high value of KDAS at 30 hours. Although the peaks of the curves for the intermediate and high values appear at about the same y-position, the carrying capacity of the wave for the intermediate value is clearly greatest. SSJS‘Var In vahd.Th meES coefficien values, wf l 1.32x10 Concenu Figure 3 source Concent 10 and i 87 5.5.5 Variation of chemoattractant concentration In this group of simulations, the initial conditions for the chemoattractant S were varied. The base case set of initial conditions for all of the chemical compounds is given in Table 5. For all simulations, the east and west reservoirs were sealed (the mass transfer coefficient was set to 0). The parameters for Population A were set to the base case values, with the exception of XOAs, which had a value of 0.05 cmZ/h. Table 5. Base case initial conditions for competition simulations. Compound and location Value (g/cm3) S in arena, south reservoir 0 S in north reservoir 1.32><10"4 Q in arena, N and S reservoirs 1.32><10'5 H in arena, N and S reservoirs 4.6x10’4 P0pulation A, Population B 5x10'6 The concentration of S in the north reservoir was varied from 1.32><10‘6 to 1.32><10'2 g/cm3. Figure 35A and Figure 35B show the dynamic competition factor for six concentrations of S in the source reservoir. Like the results shown in Figure 34A and B, Figure 35A and B also show a reversal in the \[l trend above a certain chemoattractant source concentration. For values of S S 1.32x10’6, w' is always negative. As the concentration is increased, the Sign of w‘ changes after a certain amount of time (between 10 and 20 hours) in each simulation. Above a certain concentration, however, the u! trend begins to « the trade-c Fig concentra‘ mmgnm time-pom 35C and i is that as concentr: Opposing simulatir conditio: for three indicate. growth 1 Wall of Caflying 88 begins to decrease with increasing concentration. This effect can again be explained by the trade-off between the maximum flux and saturation of receptors. Figure 35C and Figure 35D show the maximum obtainable flux for the six concentrations. Note that the y-axis has a different scale in the two figures. To avoid using simulations where the cells had encountered the wall of the chamber, the 28 hour time-point was used for all graphs that show a trend at a single time. The graphs in Figure 35C and Figure 35D have been split because the scales are so different. The overall trend is that as the concentration of S increases, so does the maximum attainable flux. Figure 35E shows the tbs curves for the four lowest concentrations. At the higher concentrations, the (1)5 values were extremely close to O at every position. The same opposing trends that were observed in the response to KDAS are apparent in this set of simulations. The level of S that minimizes saturation and maximizes flux imparts the best conditions for survival to Population A. In Figure 35F, the average concentration curves for three concentrations are shown. At the lowest concentration, a wave does not form, as indicated by the lack of a local maximum in cell concentration near the source side of the growth region. At the intermediate concentration, a strong wave has traveled nearly to the wall of the chamber. At the highest concentration, a wave is apparent, but has a low carrying capacity and has not moved as far from the inoculation point. Ar: N.Ee .39 an... x .35". 89 wt) o 10 20 30 40 Time (h) 7 C 1.32005 9 6 " .5" '2 Z: ‘3' 5 - ‘t‘ E l e o 4 o ‘5 a, i 3 v v ‘2: 3 " "b 1- 1- X 2 .. X 2 l a IL 1 T “- o . . o 1 2 a 4 5 o 1 2 3 4 x-Position (cm) x-Position (cm) 1 1.6 1.3200“ F V 1.4-4 08- E 2 l ' i x .... c .9 0.6 i 6.6x10'° ‘5 1. b 6“ 9- i\/ 3 0 0'8“ 0.4 i 1.32x10‘5 5 32 O 0.6 -- o a: 0.2 “ g o_4 ._ 1.32x10“ 5 0.2g 0 4—— r a 4. r o 1 2 3 4 5 o . . if . 0 1 2 3 4 5 X-POSitiOfl (cm) y-Position (cm) Figure 35: Variation of S source concentration. Time=28 h for all figures. 5.5.6 Var A microbial S-source Altemati‘ chemoatt (101)) in A was chemor reSpon enter 2 Compe 90 5.5.6 Variation of inoculation conditions A large variety of inoculation patterns can be envisioned that could influence the microbial competition. For example, one population could be inoculated further from the S-source than the other, possibly disadvantaging the population further from the S-source. Alternatively, two different inoculation points equidistant from the source of the chemoattractant S could be used. Figure 36. Pop. A inoculated farther from S source than pop. B. In Figure 36, a the first approach was tested. The source of S is from the north (top) in each image. The cell properties are the same as used in Section 5.5.5. Population A was inoculated further from the S—source, but directly in line with Population B. The chemotactic wave in Population A still was able to form, although it may be more in response to the chemoattractant Q than to S. As the S gradient develops over time, it will enter areas of the DGC where Population A is present, and Population A may gain the competitive advantage. i discusse moved t Populati l in Figu POpulat (approx indicati POpulat COmpet WCH-es test if 91 Figure 37. Pop. A inoculated to the side and away from S source. Figure 37 shows the results for a combination of the two inoculation classes discussed. Population B has been moved closer to the S-source, while Population A was moved to one side, away from Population B. Again the chemotactic wave is apparent, and Population A is able to move into a larger area of the chamber than Population B. Figure 38, which shows the dynamic competition factor curves for the simulations in Figure 36 and Figure 37, indicates that in both figures, the S-gradient is giving Population A some advantage, because it is only after S reaches Population A (approximately 10 hours), that the dynamic competition factor begins to increase (rp'>0), indicating that Population A is gaining mass at a faster rate than B. It is apparent from these few test situations that the inoculation pattern of the two pOpulations can play an important role in determining which population will have the competitive advantage. Another situatiOn that could be tested is to have one population well—established in the DGC, and introducing the other at a much lower concentration, to test if it can successfully "invade" and be coexistent with the initial population. The results of the time 1 5.5.7 i compo being Popula nutrier hl'potl given 92 results of all of these variation of inoculation point simulations are highly dependent upon the time point being studied. 1 0.3 -— 0-6 " Figure 36 A ,o 4-0 ,1 V _ 0.4 -- - 0.2-- t Figure 37 0 l l i r r o 5 10 15 20 25 30 Time (h) Figure 38. Dynamic competition factors for two inoculation patterns. 5.5.7 Variation of nutrient initial conditions Another hypothesis that was explored with the model was that if the nutrient compound was diffusing in from one reservoir with the chemoattractant, as opposed to being uniformly distributed through the gel, then the chemotactic response of a pOpulation might give it the competitive advantage by moving that population toward the nutrient, while the other was left behind in a low nutrient environment. To test this hypothesis, a simulation was performed where the two populations had the parameters given in Table 4, except XOAs=O.05 cmZ/h, VAH=O-5 hJa and VBH=O-7 114' The nutrient was present i was pre: same as simulati 39C sho 93 present in the north reservoir (x = 5 cm) at a concentration of 8.6><10'4 g/cm3, and no H was present in the gel initially. The chemoattractant and cell initial conditions were the same as given in Table 5. Figure 39A shows the dynamic competition factor for this simulation, and Figure 39B shows the average concentration curves at time = 30 h. Figure 39C shows the centerline profiles (y=2.5 cm) of the H and S compounds. 3 i-i )< 1 O3 gg/ksrvw HCVer Populz 393. ti 0 O 94 3.5 6‘ 0.8 E PopulationA 3 T a 0‘7 T ----- PopulatlonB L “i: 0.5-i 2.5 — F x L g 0.5 «i A 2- 3 =3 S 0.4-- "‘ — a 15 i 8 0.3i 8 1 0 o 0.21. 8’ 0.5 a I- 0.1 . l s 1 ‘ < o + + : 4| + 0 ~41; . 0 5 10 15 20 25 30 0 1 2 3 4 5 Time (h) Position (cm) 1 0.12 0.8 L - - H profile " S profile _ 009 (V) (‘0 E 8 0.6 T 5 CD C) no _ 0.06 ”O 1- 1— x 0.4 "L' x I U) _ 0.03 0.2 i— I” 0 d I. _____*--‘---"T- O Position (cm) Figure 39. Competition for nutrient diffusing in at x = 5 cm. Figure 39A illustrates that even though Population B has a growth advantage, it never obtains a higher mass than Population A. This is explained by the fact that POpulation B never experiences a high concentration of the nutrient. As shown in Figure 39B, the chemotactic response of Population A has allowed it to move toward the source of nutrie complete mass of likely is consume 1 pronoun was a f: exist. cl examplt discuss 5.6 F‘ methor should single Cheme and th Studie hOWe‘ 95 of nutrient and chemoattractant. Figure 39C illustrates . that the nutrient has been completely consumed by P0pulation A, and is not available to Population B. The small mass of Population B that appears close to the source of S and H in Figure 39B most likely is the product of a small amount of growth that occurred before Population A fully consumed the H that had diffused to the center of the chamber. In this simulation, the chemotactic response of Population A gave it an even more pronounced advantage in that u: never had a value less than 1, even though Population B was a faster grower. In a system where gradients of both chemoattractants and nutrients exist, chemotaxis may be extremely important in determining the competitive winner. An example of a real system where these conditions might occur is bioremediation, as discussed in Section 7. 5.6 Future work on competition The next step needed in the competition studies is to develop experimental methods to corroborate the modeling predictions. The focus of the experimental work should be on finding mutant bacteria that would differ from the main strain only in a single transport property, such as random motility, or response to a single chemoattractant. Many transport mutants of E. coli have been previously isolated and cataloged, and the genes involved in the specific mutations are well defined. Emerson et al. (1994) studied the growth of such E. coli mutants in the DGC. At the time of this writing, however, Pseudomonas KC (see Section 7) was the focus of this portion of the research project. l efforts tc working mutants desired 1 possess mutant 2 after the marker competf Work) 1 Escherr more ( mediur asparta grown. Optical l:igure 96 project. Well-characterized mutants of Pseudomonas KC are not currently available, but efforts to develop them are underway in Dr. Craig Criddle's laboratory. Petty Setiawan, working on a project that was part of this thesis work, screened 120 Pseudomonas KC mutants for their motith and chemotaxis parameters. None were found that have the desired prOperties that would make them useful in the competition experiments. One way to perform these experiments would be to select transport mutants which possess a selection marker, such as antibiotic resistance or fluorescence. A mixture of the mutant and non-mutant strains could then be inoculated in a swarm plate or the DGC, and after the growth and motility patterns had formed, samples could be taken. The selection marker would allow each strain to be counted individually, and the outcome of the competition could be quantified. A preliminary competition experiment has been performed (Petty Setiawan's work) to begin to study the effects of chemotaxis on competitions. In this experiment, Escherichia coli was competed against Pseudomona stutzeri strain KC (see Section 7 for more details on Pseudomonas KC). Swarm plates were poured that contained M9 medium at pH 7.5. The plates contained 2 mM glycerol as the nutrient source and 0.1 mM aSpartate as the chemoattractant. Two plates were inoculated with 20 ill of either actively growing Pseudomonas KC or E. coli culture that had been adjusted to have the same Optical densities. The plates were incubated at 27 °C for 48 hours. These plates, shown in Figure 40, allowed the response of the individual populations to be observed. Fi density for 48 sample diluted Were i1 KC C0 97 Mixture of two populations Figure 40. Competition between Pseudomonas KC and E. coli in swarm plates. In a third swarm plate, a 20 ul aliquot of an equal mixture, based on optical density, of the two populations was inoculated in a swarm plate and incubated at 27 °C for 48 hours. An image of this swarm plate is also given in Figure 40. After 48 hours, samples were taken at 5 locations in the competition swarm plate. The samples were diluted by a factor of l:10,000 and plated on nutrient agar plates. The nutrient agar plates were incubated for two days at 27 °C. Individual colonies of E. coli and Pseudomonas KC could be differentiated by their unique morphologies, and counted. The results are shown it plate. results The P5 the €te illiO arr gradiei (glycei Was in 98 shown in Figure 41. The origin of the x-axis corresponds to the center point of the swarm plate. onn Concentration x 104 (cells/ml) Position (cm) Figure 41. Colony counts in competition experiment. The E. coli remained primarily in the center of the swarm plate. Although the results are not symmetrical, the Pseudomonas KC have moved outward from the center. The Pseudomonas KC may be at a competitive advantage in this case if the resources in the center of the plate are fully consumed. The Pseudomonas KC would be able to move into areas at the edge of the growth front where resources are still plentiful. More work needs to be done with experiments such as this one. For instance, gradient measurements would be useful in determining if resources such as the nutrient (glycerol) were still available in the center of the plate. Also, in this experiment, no effort was made to control the relative growth rates of the two populations. 6.5fldl a M. S. were pe Mutaget Chambt Bioclzei DAHP metabo feedbar tested 1 had in anaIOg away f lid Wa: for 10 that w The cl glucos 6. SELECTION OF MUTANTS 6.1 Experimental work This portion of the research was completed as a collaboration with Mark Mikola, a M. S. candidate in Dr. Worden's research group. The experimental portions of the work were performed entirely by Mr. Mikola, and are detailed in Mikola (1996) and in "In-situ Mutagenesis and Chemotactic Selection of Microorganisms in a Diffusion Gradient Chamber" by M. R. Mikola, M. T. Widman, and R. M. Worden, submitted to Applied Biochemistry and Biotechnology. To summarize, it was desired to obtain a strain of E. coli whose isozyme of DAHP synthase, designated AroF, was not inhibited by L-tyrosine, a product in the metabolic pathway. To achieve this goal, a population of E. coli that experienced feedback resistance from tyrosine was inoculated into a DGC. The E. coli had previously tested positive for chemotaxis toward glucose in a separate DGC experiment. The DGC had in its source flask 125 ttM m-fluorotyrosine (m-FT), a non-metabolizable tyrosine analogue, and 5 mM glucose. The bacteria only grew in the portion of the DGC farthest away from the source of m-FT (see Figure 42A). After the initial growth period, the DGC lid was removed, and the DGC was exposed to ultraviolet radiation at a distance of 60 cm for 10 seconds. The lid was replaced and the incubation continued. Blooms of mutants that were less inhibited by the m-FT appeared after the mutation process (Figure 42B). The chemotactic response to the glucose gradient drew the mutants toward the source of glucose and m-FT, as seen in Figure 42C. 99 -— Fig 6.2 modifie cell ba' Also. a by but] Where consu satur; 100 Figure 42. Mutant selection by chemotactic response (source at top of pictures). 6.2 Modeling selection of mutants The mathematical model of the DGC system presented in Section 4.2 was modified and used to simulate strain selection in the DGC. For this application, a second cell balance was added to account for both the original (A) and the mutant (B) strains. Also, a balance was added to account for diffusion and cellular uptake of the inhibitor (P) by both strains: :92_ 8t _ vaFP vbPP ( 46) a ”[7 Cap + P C“. + P DPVZP— where Dp is the diffusion coefficient of the inhibitor; Vap and pr are the specific consumption coefficients for Populations A and B consuming P; and Cap and Cbp are the saturation constants for consumption of P by Populations A and B, respectively. chemoat growth 1 terms cl Ollis, 15 where respont sensitii pOpula by Where uninh motili 101 The chemotactic sensitivity coefficients for both populations responding to the chemoattractants S and Q, the random motility coefficients, and the maximum specific growth rates were modified to incorporate the inhibition effects. The form of the modified terms chosen is analogous to noncompetitive inhibition in enzyme kinetics (Bailey and Ollis, 1986). The modified chemotactic sensitivity is given by , _ x... (47) where x'Oij is the inhibited chemotactic sensitivity for the ith population (i=a or b) responding to the jth chemoattractant (j=S or Q), XOij is the uninhibited chemotactic sensitivity, and Kfijx is the inhibition constant for the chemotactic sensitivity of the ith population to the th chemoattractant. The modified random motility coefficient is given by I “i ( 48) where u'i is the inhibited random motility COCfflClCI'lt for the 1t population, it, is the uninhibited random motility coefficient, and K1,l1 is the inhibition constant for the random motility of the ith population. The inhibited maximum specific growth rate is given by where v H, till 1 specific cell bal where includt Where the yi (Widi Equat trend: funct Calcu I vii-l ( 49) where v'iH is the inhibited maximum specific growth rate for the ith population growing on H, win is the maximum specific growth rate, and KEV is the inhibition constant for the specific growth rate of the ith population. In terms of the variables defined above, the two cell balance equations are given by ( 50) where u, is the ith cell population (either a or b). The nutrient balance was also modified to include the inhibited maximum specific growth rates: 8H ———=D VZH— at ” _aiaiga._____véaH _u_,_ (51) C.H+H Yay CbH+H YbH where DH is the diffusion coefficient for H, Cm is the half-saturation constant, and YiH is the yield coefficient. Values for the modeling constants were taken from the literature (Widman et al., 1997) insofar as possible. Values for the inhibition constants included in Equations ( 47)-( 49) were chosen that gave reasonable agreement with the experimental trends. Figure 43 shows a time sequence of contour plots depicting cell concentration as a function of position. The contour lines show constant concentration isoclines, as calculated by a built-in Matlab program, The chemoattractant (glucose) gradient is indkxne Thereis thetop observe 103 indicated by gray shading, with lighter gray indicating a higher glucose concentration. There is no m-FT in this simulation. The source reservoir is on the side correSponding to the top of each figure. Bias of the cell’s migration in this direction is comparable to that observed experimentally. 104 Increasing Time Figure 43. E. coli responding to chemoattractant gradient. shown evohni (nuytt unhcat 43.is inhhnt DCKl paant Popul bound The - cheni secon the g reser and . inoct POpU 105 The ability of the mathematical model to reproduce the experimental trends shown in Figure 42 was evaluated. Figure 44 shows the predicted growth-pattern evolution in which gradients of both a chemoattractant and an inhibitor are applied, and only the population sensitive to the inhibitor is present. In this figure, the gray shading indicates the inhibitor gradient. The effect of chemotaxis, which was significant in Figure 43, is overwhelmed by the effect of the inhibitor. As the concentration gradient of the inhibitor becomes established, cell growth occurs predominantly in the sink end of the DGC. This trend was observed experimentally in Figure 42. The base case set of parameters for all of the simulations shown is given in Table 6. In these simulations, Population A was the inhibited population, while Population B was the mutant. The boundary and initial conditions are similar to those presented in Sections 4.2.1 and 4.2.2. The inhibitor concentration (P) in the source reservoir was 0.0132 g/cm3. The chemoattractant concentration (S) in the source reservoir was 0.000132 g/cm3. The second chemoattractant concentration (Q) in the source and sink reservoirs and initially in the gel was 0.0000132 g/cm3. The nutrient concentration (H) in the source and sink reservoirs and initially in the gel was 0.000046 g/cm3. The value of uao was 3><10‘6 g/cm3, and ubO was 3><10'7 g/cm3. The gradients were allowed to initiate for 12 h before inoculation of Population A. Population B appeared in the simulation 15 h after Population A. The simulation continued for 25 h after the appearance of Population B. The shown in I evolution in only the p0j indicates [iii 43, is oven inhibitor be DGC. This parameters Population boundary; The inhib chemoattri second Chi the gel w reservoirs and ub0 inoculatic POPUlatio 105 The ability of the mathematical model to reproduce the experimental trends shown in Figure 42 was evaluated. Figure 44 shows the predicted growth—pattern evolution in which gradients of both a chemoattractant and an inhibitor are applied, and only the population sensitive to the inhibitor is present. In this figure, the gray shading indicates the inhibitor gradient. The effect of chemotaxis, which was significant in Figure 43, is overwhelmed by the effect of the inhibitor. As the concentration gradient of the inhibitor becomes established, cell growth occurs predominantly in the sink end of the DGC. This trend was observed experimentally in Figure 42. The base case set of parameters for all of the simulations shown is given in Table 6. In these simulations, Population A was the inhibited population, while Population B was the mutant. The boundary and initial conditions are similar to those presented in Sections 4.2.1 and 4.2.2. The inhibitor concentration (P) in the source reservoir was 0.0132 g/cm3. The chemoattractant concentration (8) in the source reservoir was 0.000132 g/cm3. The second chemoattractant concentration (Q) in the source and sink reservoirs and initially in the gel was 0.0000132 g/cm3. The nutrient concentration (H) in the source and sink reservoirs and initially in the gel was 0.000046 g/cm3. The value of uao was 3><10'6 g/cm3, and ubO was 3x107 g/crn3. The gradients were allowed to initiate for 12 h before inoculation of Population A. Population B appeared in the simulation 15 h after Population A. The simulation continued for 25 h after the appearance of Population B. 106 Table 6. Parameter values for inhibition model Parameter Value (i=a or b) p, 0.010 cm2 h'l XOiQ 0.08 chh" m, 0.08 cmth CiH 4.08x10'6 gH cm'3 i CiQ 6.70x10'8 gQ crn‘3 CS 5.50x10'5 gs crn‘3 DH 0.01 or? hj DQ 0.033 crthr DS 0.033 cmfh’l Dp 0.033 cmfh'l KDiQ 3.30x10'5EQ crn'3 KDis 2.00x10'fgs crn'3 v,“ 0.35 h~I ViQ 0-02 gQ gr:1 h-1 vis 0.60 gs g,,'1 h’1 YiH 0‘50 gu/gH Klajx 1.0x10'6 g cm‘T K1,, 10x10T g crnT Km 1.0x10'6 g cm: Figui 107 Increasing Time A Figure 44. E. coli responding to inhibitor gradient and chemoattractant gradient. Figt B) that is it the inhibitr significant 3 lower cor effect of th moves prt predomina experimen Th1 favorable that lose I chemoattr: method sl inhibtion i the trends be elucida 108 Figure 45 is a simulation of a mutation event giving rise to a mutant (Population B) that is insensitive to the inhibitor. As in Figure 44, Population A, which is sensitive to the inhibitor, grows preferentially near the sink. Population B takes longer to appear in significant concentration, because this population was initiated after Population A, and at a lower concentration than Population A, to simulate the UV mutagenesis. The beneficial effect of the chemoattractant in separating the two populations is evident, as Population B moves preferentially toward the source reservoir, while Population A remains predominantly near the sink reservoir. These trends are similar to those observed experimentally in Figure 42. The model allowed the effectiveness of the method to be explored under less favorable conditions than were experienced in the experimental work, such as mutants that lose their chemotactic re3ponse, or that appear in areas of the DGC far from the chemoattractant source. The model predicts that even in these worse case scenarios, the method should still allow the desired mutants to be selected. The noncompetitive inhibtion form used in the modeling (Equation ( 47), for example) adequately reproduced the trends of the experiment, but other forms may be more physically realistic, and may be elucidated through further experiments and modeling. 109 Increasing Time \\.=. A Figure 45. Simulation of experiment shown in Figure 42. The r asimulation of Populatio two populat predominan to Q. T0 es repeated wf chemotaxis potent the t 1 10 The model can also be used to test other situations. For example, Figure 46 shows a simulation of the case where the mutation occurred in the middle of the growth pattern of Population A, and the mutant (Population B) was not chemotactic to S. In this case, the two populations were still able to be separated, but in this case the separation is based predominantly on the influence of the inhibitor gradient and spreading due to chemotaxis " to Q. To estimate the benefit obtained by chemotaxis toward S, this simulation was then repeated with chemotaxis to S reinstated. The results, shown in Figure 47, indicate that chemotaxis toward S does further enhance the rate of separation. Presumably, the more potent the chemoattractant (i.e., the higher the x'mj value) the greater the enhancement. 111 Increasing Time Figure 46. Mutant appears far back in population. 112 Increasing Time Figure 47. Mutant now chemotactic to S. 7. BIORI The environmen e101,. 1993: references, or consider might sugg that severa through so: Rec porous me Walk mod the ratio ( porous mt column ti. both the r replaced eCOiOgica roots by j and Chen- 7. BIOREMEDIATION The transport of bacteria through porous media, specifically soil or aquifer environments, has been studied in both experimental and modeling systems (Abu~Ash0ur, et al., 1993; Sarkar, et al., 1994a; Sarkar et al., 1994b, for example). In a number of these references, the chemotactic movement of the bacteria has been ignored (Tan, et al., 1994) or considered too complex to include, even though results have been recorded which might suggest chemotaxis is important. For example, Abu-Ashour et al. (1993) reported that several field and laboratory experiments found that the average velocity of bacteria through soil was faster than the velocity of tracers and of the groundwater. Recently, a few studies have addressed the role motility plays in transport through porous media. Duffy et al. (1995) attempted to adapt the homogeneous media random walk model to a porous media. They utilized a tortuosity factor that was proportional to the ratio of the random motility in bulk liquid to the effective random motility in the porous media. Their model gave good qualitative agreement with an experimental sand column through which Pseudomonas putida swam. Barton and Ford (1997) showed that both the random motility coefficient and the chemotactic sensitivity coefficient could. be replaced with effective values that incorporate the effect of the porous media. An ecological system in which chemotaxis through soil is important is the colonization of roots by Rhizobium meliloti. Soby and Bergman (1983) demonstrated that active motility and chemotaxis were necessary for efficient spreading of the bacteria through soil. 113 The: an in situ phenanthren the areas of They blame sands. reter contaminar presented 2 model inc experimen PSt tetrachlori chloroforr other mic microbe ' €Xperime 1996). 1 bioaugm. Schoolcr as an elt aCccptor COmami: 114 The effective transport of bacteria through the soil is an important requirement for an in situ bioremediation project. Devare and Alexander (1995) found that a phenanthrene-metabolizing Pseudomonas sp. could adequately remove phenanthrene in the areas of a soil column in which the bacteria was inoculated, but not in other areas. They blamed this result on bacterial retention by the clay soil, and noted that in aquifer sands, retention was much less. Their conclusion was that transport of bacteria to the contaminant was very important to the success of bioremediation. Bosma et al. (1988) presented a model for bacteria consuming xenobiotic chemicals in a soil column. Their model incorporated chemotactic movement, and its results compared well with experimental data. Pseudomonas stutzeri KC is a denitrifying bacterium able to degrade carbon tetrachloride (CT) into carbon dioxide and nonvolatile products without the production of chloroform (Criddle et al., 1990; Dybas et al., 1995; Mayotte et al., 1996). In contrast, other microbes typically convert CT to chloroform. The outstanding potential of this microbe for bioremediation of CT spills has been demonstrated both in shaker-flask experiments (Tatara et al., 1993), and a model aquifer system (Witt 1994; Mayotte et al., 1996). The Michigan Department of Natural Resources is sponsoring a major bioaugmentation field experiment with this organism in a CT-contaminated aquifer near Schoolcraft, Michigan. In this experiment, acetate is periodically injected into the ground as an electron donor for the reaction, and the naturally occurring nitrate is the electron acceptor. The Pseudomonas KC is introduced into the aquifer soil, and the flow of contaminated water is directed through a "biofence" of the bacteria. Chet in the labo packed witl equal to tha indicated ti front was n 15 cm/day. 7.1 Expt 7.1.1 Swa Stu KC t0 var were prep; adjusted tr SOlutionsr into petri Pseudmm exlitifimer SCientific Fi levels of: nitrate c0 1 15 Chemotactic movement of Pseudomonas KC has been observed in sand columns in the laboratory (personal communication with Mike Witt). The sand column was packed with aquifer sand, and supplemented with acetate. The nitrate concentration was equal to that observed in the Schoolcraft aquifer. Preliminary experimental measurements indicated that the chemotactic movement was in response to nitrate gradients. The cell front was moving at a velocity approximately 8 cm/day faster than the water velocity of 15 cm/day. 7 .1 Experimental results 7 .1.1 Swarm plates Studies were initiated to characterize the chemotactic response of Pseudomonas KC to various possible chemoattractants. Solutions of Medium D (Tatara, et al., 1993) were prepared with varying levels of acetate and/or nitrate. The pH of the medium was adjusted to 8.2. Swarm plates (Petty Setiawan's work) were made using the Medium D solutions supplemented with 0.25% high strength agar gel. The gel mixture was poured into petri dishes and allowed to solidify. A 20 til aliquot of actively growing Pseudomonas KC culture was inoculated in the center of the plate. For the anaerobic experiments, the plates were incubated in a GasPak 150TM Anaerobic System (VWR Scientific). For aerobic experiments, the plates were incubated on the lab bench. Figure 48 shows the chemotactic ring patterns formed in response to varying levels of nitrate, at a constant acetate concentration of 1 g/l. In the image on the left, the nitrate concentration was 1 g/l, and on the right, the nitrate concentration was 2.5 g/l. A11 swarm plate to l g nitre Pseudomom acetate to 11 image of F occurred. C inner ringt consumptit may be du- experirr the left 1 l6 swarm plate images were taken 48 hours after inoculation. The mass ratio of l g acetate to 1 g nitrate is the balanced ratio for consumption of the two chemical species by Pseudomonas KC (personal communication with Dr. Craig Criddle). When the ratio of acetate to nitrate was equal to one, only a single ring appeared, as shown in the left-hand image of Figure 48. In the right-hand image, when nitrate is in excess, a double ring occurred. One hypothesis is that the outer ring is due to consumption of acetate, and the inner ring due to the use of the remaining nitrate, or an intermediate such as nitrite, for the consumption of another carbon source, such as endogenous metabolism. The outer ring may be due to acetate, nitrate, or a combination of both. Figure 48. Chemotactic rings of Pseudomonas KC at varying nitrate levels. In the next experiments, the acetate concentration was varied. Figure 49 shows experiments at a constant nitrate level of 50 mg/l, and acetate varying from 100 mg/l on the left to 1000 mg/l on the right. In the low acetate image, there is no noticeable movement i miniml let considerablt Figure 48. I hand image image 100( the many ' appeared v bubbles m denitriftca‘ 117 movement from the inoculation condition. Both growth and motility seem to be at a minimal level. In the higher acetate image, however, the high density ring indicates considerable growth, but the diameter of the ring indicates that motility is slower than in Figure 48. In the next set of experiments, the nitrate level was set to 500 mg/l. In the left- hand image in Figure 50, the acetate concentration was 100 mg/l, and in the right-hand image 1000 mg/l. At 100 mg/l acetate, only a single ring forms. In the 1000 mg/l image, the many bright white dots are bubbles that appeared in the gel. The bubbles only appeared when the acetate to nitrate ratio was greater than one. We hypothesize that the bubbles may contain nitrogen that forms when excess acetate is present to complete the denitrification reaction shown in Equation ( 52) Figure 49. Nitrate concentration = 50 mg/l, varying acetate. Pseudomc acetate is acceptor t nitrogen, mEdium. found ti r2iiiging fOrrned. at the h] Pseudo 118 Figure 50. Nitrate concentration = 500 mg/l, varying acetate. Pseudomonas KC is known to preferentially use nitrate before consuming nitrite. When acetate is in excess, however, the Pseudomonas may utilize the nitrite as an electron acceptor to further consume the remaining acetate. This would result in the formation of nitrogen, and possibly the bubbles. No; —> No; —> N2 (52) In aerobic experiments with the Pseudomonas KC, no nitrate was added to the medium. Instead, oxygen was used as the electron acceptor. In results not shown, it was found that no ring formed at an acetate concentration of 100 mg/l. For concentrations ranging between 1000 mg/l and 1750 mg/l, a ring of approximately 4.7 cm in diameter formed. At a concentration of 2500 mg/l, no ring was observed. These results suggest that at the higher concentrations of acetate, the chemotactic receptors may be saturated. Additional experiments were performed to test the aerobic response of Pseudomonas KC to attractants other than acetate. Medium M9 (Maniatis et al., 1982), adjusted to medium we more cells experiment supplemen carbon sou 119 adjusted to pH 7 .5, was used for the medium. In the left-hand image in Figure 52, the M9 medium was supplemented with 1 mM glucose. A distinct chemotactic ring formed, but more cells seemed to remain in the center of the plate than in the acetate/nitrate experiments shown above. In the right-hand image of Figure 52, the medium was supplemented with 0.1 mM aspartate as the chemoattractant, and 2 mM glycerol as the carbon source. In this instance, a fairly wide, diffuse ring was observed. Figure 51. Pseudomonas KC responding to glucose (left) and aspartate(right). 7.1.2 DGC experiments Experiments were performed in the DGC to test the response of the Pseudomonas KC to applied gradients of acetate and nitrate. Medium D-, containing 3.5 g/l K2HP04, 1.24 g/l KH2P04, and 1.0 g/l (NH4)2SO4 was used as the minimal medium for the experiment. The pH of the medium was adjusted to 8.1 with KOH. After autoclaving, 2 mid of l M to test the g experiments comparable contained 1 minimize t chemotacti degradatio: An nitrate ant addition tt sink flask and 0.25‘. C011 grow toward ti nearest t alilieared Source 1 Samples the DG( 151118 120 m1/1 of 1 M MgSO4 was added to the medium. Shake flask experiments were performed to test the growth characteristics of Pseudomonas KC in this medium. The shake flask experiments showed that the maximum growth, as measured by optical density, was comparable to growth in the complete Medium D recipe. The Medium D- recipe contained fewer components, which was useful not only for simplicity, but also to minimize the number of consumable chemical components that could possibly elicit a chemotactic response. The effects of using Medium D- on carbon tetrachloride degradation were not studied as part of this work. An initial DGC experiment was performed to test the chemotactic response to nitrate and acetate in Medium D- (results not shown). The source flask contained, in addition to 800 ml of Medium D—, 5 mM sodium acetate and 5 mM sodium nitrate. The sink flask contained 800 ml of Medium D—. The arena of the DGC contained Medium D- and 0.25% high-strength agar gel. After 28 hours of growth, a pattern similar to the E. coli growth patterns shown in Figure 13, had developed, with a definite bias in movement toward the acetate/nitrate source. A 0.5 ml sample was removed from the gel at the point nearest to the source reservoir. A similar sample was removed from the gel at what appeared to be the center of the chemotactic wave moving toward the acetate/nitrate source reservoir. The samples were injected into shake flasks and incubated. These samples were taken in an effort to isolate bacteria that were optimized for chemotaxis in the DGC environment. A second DGC experiment was performed with the same parameters. This time, a 15 til aliquot of the cell culture taken from the spot nearest the wall in the previous experiment 22 hours b: Figure 52. inoculation 121 experiment was used to inoculate the chamber. The gradient was allowed to initialize for 22 hours before inoculating the bacteria. The results of this experiment are shown in Figure 52. The times shown in the corners of the images correspond to the time after inoculation of the bacteria. The acetate/nitrate source is at the bottom in each image. 122 Figure 52. Pseudomonas KC responding to acetate and nitrate in a DGC. 7.1.3 Mo‘ A otherwise objects. Tr experimet sodium at KOH. Ti connector solidify, Figure 53 across th ends oft the wave 123 7 .1.3 Movement of wave around obstacles A first logical step in understanding the effect of solids on chemotaxis in an otherwise homogenous medium would be to study a medium which has one or a few objects. To test the response of the chemotactic wave to obstacles in its path, swarm plate experiments were performed (Laura Booms' work). Medium D- supplemented with 5 g/L sodium acetate, 3 g/L sodium nitrate, and 0.25% agar was adjusted to a pH of 8.2 with KOH. The plates were poured, and before the agar solidified, two sterilized tubing connectors with sealed ends objects were placed in the gel. The gel was allowed to solidify, and then 20 111 of Pseudomonas KC culture were inoculated in the center. In Figure 53A, taken 26 hours after inoculation, the wave has moved approximately halfway across the tubing connectors. In Figure 53B (42 hours), the wave has moved beyond the ends of the connectors. The objects did not appear to perceptibly inhibit the movement of the wave. Figure 53. Movement of wave around obstacles in swarm plates. To performed diameter 1 the path 0 P5 and the st in a micrc 15 pl 0ft the gradit A The time the botto in Figure experimt will not 124 To test how a large object would affect the chemotactic wave, an experiment was performed under the following conditions in a DGC. A cut-off glass test tube 1 cm in diameter was inserted into the gel, close to the source reservoir, to serve as an object in the path of the chemotactic wave. Pseudomonas KC were grown for 3 days in a 100 ml flask containing Medium D- and the same ratio of acetate and nitrate as in the source flask. 6 m1 of culture were spun in a microcentrifuge, and the resulting pellets were resuspended in 0.5 ml of fresh media. 15 1.11 of the resuspended culture were inoculated into the center of the arena 5 hours after the gradient was initiated. A time-sequence of image captures from the experiment is shown in Figure 54. The time increases from left to right, and t0p to bottom. The source reservoir is located at the bottom in each picture. The edge of the wave just encountered the tube in Figure 54A. In Figure 54D the cells have passed the object and have rejoined on the other side. This experiment indicates that, at least in this very simple single-object situation, an obstacle will not catastrophically disrupt the chemotactic movement of the cell population. 125 Figure 54. Chemotaxis of Pseudomonas KC around an object. 7.2 Fut Ft hypothesi focus on a wall of DGC wa cm from into the allowed solidifie Experirr designer glass be Structur could b. 7.2 Future work on bioremediation Further experimental and modeling work will be necessary to prove the hypothesis that chemotaxis can enhance a bioremediation effort. First, studies should focus on the effects of a porous medium on the chemotactic response. A method to create a wall of glass beads, sand, or soil in the DGC has been developed (see Figure 55). A DGC was modified to allow two specially cut microscope slides to be placed 1 cm and 2 cm from one wall. Glass beads were placed between the slides, and the gel was poured into the chamber. Small gaps between the slides and the bottom cover of the DGC allowed the gel to penetrate under the slides and up into the beads. After the gel solidified, the slides were carefully removed, and the three gel regions combined. Experiments to test how the bacteria are transported through the porous wall should be designed and performed with and without chemoattractant gradients present. Initially, glass beads of known diameter could be used to give a system with well—defined pore structure. Next, a more environmentally realistic wall material, such as aquifer sands, could be used. 126 v Chemo source contami had a it used to model 1 model Worder pOI‘OUS Carbon 127 Microscope slides Chemoattractant F3611 1 ' source Inocu atron pornt Figure 55. Set-up to make glass-bead wall in DGC. A version of the mathematical model was developed that included a balance for a contaminant. The balance was similar to the nutrient and chemoattractant balances, but had a term that accounted for disappearance of the contaminant. This model could be used to predict the enhanced contaminant degradation when chemotaxis occurs. The model was never fully tested and debugged, and suitable parameters for the contaminant model were not found. A new cellular dynamics model has been proposed by Drs. Worden and Lastoskie which would allow for predictions of the behavior of cells in a porous medium. For the Pseudomonas KC system specifically, a method to quantify the amount of Carbon tetrachloride in a DGC will need to be developed. A control experiment with no applied c? degradatic movemen' compare 1 growth st independe Or bioremed cells witl injection possible uncontan pulse che how the (Section maximiz 0Plimize zones. 1 bioreme money ( 128 applied chemoattractant gradient should be performed to measure the amount of degradation of carbon tetrachloride that takes place in the absence of chemotactic movement. Then an experiment with an applied chemoattractant gradient should be run to compare the degradation rate to the control experiment. Factors such as the amount of growth should be accounted for in the analysis, so that the effect of chemotaxis can be independently evaluated. Other experiments to discover ways in which chemotaxis could improve a bioremediation effort could focus on the concept of controlling the dissemination of the cells with chemoattractant gradients. One hypothesis is that by carefully placing the injection wells through which the chemoattractant is introduced into the soil, it may be possible to direct the spreading of the bacteria so as to minimize the losses into uncontaminated zones. A relatively simple experiment to test this hypothesis would be to pulse chemoattractant into different reservoirs of the DGC at timed intervals, and observe how the cell population responded. In the section of this thesis on microbial competition (Section 5), the model predicted that optimal chemoattractant concentrations exist that maximize the chemotactic wave properties, such as the carrying capacity, and could optimize the delivery protocols to obtain the best transport of cells to the contaminated zones. The mathematical model should prove to be an extremely useful tool to explore bioremediation, and, through modeling simulations, considerable experimental time and money can be saved. 8. SUM Tl engineeri response. nricrobio diffusion chemoat' mathem: having I predictit involve: respons mathen the pre advant technic compe eXploi could inhibi 8. SUMMARY The focus of this research was on microbial chemotaxis, and how it could benefit engineering applications. Tools were developed to aid in the study of the chemotactic response. Experimental tools included the diffusion gradient chamber, microsensors and microbiosensors, and the laser diffusion capillary assay. A mathematical model of the diffusion gradient chamber was also developed that allowed two cell balances, two chemoattractant compounds, and a nutrient compound to be modeled simultaneously. The mathematical model allowed predictions of chemotactic responses to be made without having to run time consuming experiments. The model was validated by comparing its predictions to experimental data from DGC experiments. The validation experiments involved Escherichia coli responding to gradients of aspartate. The applications that were identified that could benefit from the chemotactic response were microbial competition, selection of mutants, and bioremediation. The mathematical model was used to explore the effects of several modeling parameters on the predicted competitive outcomes. Chemotaxis was shown to give a competitive advantage to a bacterial population in certain cases in non-mixed environments. Several techniques to analyze the modeling predictions were developed, including the dynamic competition factor. For the selection of mutants t0pic, the chemotactic response was exploited to draw mutants of a bacterial species into an area in which the parent strain could not grow well. Again, the mathematical model was used to more fully explore the inhibition and chemotaxis interaction. The model predicted that the selection method 129 should wr bioremed mmdm obstacles hmMm were per benveen asanaqt 130 should work well in various scenarios. Finally, the response of an important bacterium for bioremediation, Pseudomonas stutzeri strain KC, was studied. Pseudomonas KC was found to be chemotactic to nitrate and acetate. The reaction of the chemotactic wave to obstacles also studied as a first step in examining chemotaxis through porous media. Preliminary competition experiments between Pseudomonas KC and Escherichia coli were performed. Further experiments could lead to a better understanding of the interplay between chemotaxis and competition in environments important to bioremediation, such as an aquifer. APPENDICES APPENDIX A APPEl‘ Miscella 1‘ section 1 chemotz organis: require: aliquot 56A, g gel. Th be attr oxyger and 1 define aspan multi} how the of after resea: APPENDIX A Miscellaneous chemotaxis experiments Appendix A includes selected experimental results that do not directly fit into a section of the thesis, but that are interesting or informative in the context of studying the chemotactic response. , Figure 56 shows three swarm plate experiments. In each of these experiments, the organism used was E. Coli HCB 33. M63 minimal medium, supplemented with the required amino acids and streptomycin (see Section 4.2.3), was the base medium. A 10 111 aliquot of actively growing culture was inoculated into the center of each plate. In Figure 56A, glucose was initially present at a constant concentration of 3 mM throughout the gel. The image was taken 20 hours after inoculation. The multiple waves or bands could be attributed to chemotaxis to other chemoattractants, one of which would likely be oxygen. Figure 56B, also 20 hours after inoculation, shows the response to 1 mM glucose and 1 mM aspartate. In this case, a very uniform single wave formed, that is sharply defined on the outer edge. Since at least two known chemoattractants (glucose and aspartate) were present, and possibly a third with oxygen, it might be assumed that multiple bands would occur. This did not occur experimentally in this particular case, however. It is possible that the consumption of the two chemoattractants was similar, and the chemotactic waves overlapped closely. The image in Figure 56C was taken 42 hours after inoculation. This swarm plate contained 1 mM glucose and 3 mM glycerol. Some research has reported that glycerol can be an inhibitor of chemotaxis (see Zhulin et al., 131 1997 f0 conditio Patten format 1995; Symmt 0r S. 1 substa and g Strept. inocu M63. 132 1997 for a review), which might explain the slower swarming under the glycerol conditions. Figure 56. Swarm plates with various chemoattractants. Pattern formation An interesting phenomenon that is thought to be related to chemotaxis is the formation of patterns in bacterial populations (Agladze et al., 1993; Woodward et al., 1995; Budrene and Berg, 1995). In these studies, swarm plates were used to study the symmetrical patterns, sometimes described as looking like snowflakes, formed by E. coli or S. typhimurium. The current hypothesis is that the bacteria exude a chemoattractant substance, which causes conglomeration of bacteria in patterns. A DGC experiment was set-up to test the response of E. coli HCB 33 to aspartate and glucose. M63 minimal medium was used with the required amino acids and streptomycin. The source concentrations were 3 mM aspartate and 5 mM glucose. The inoculation culture was grown in a shake flask with 10 mM glucose and supplemented M63. The gradients were allowed to initiate for 5 1/2 hours prior to inoculation. A 10 ml aliquot the cent until 1t DGC, concen' hours) produ HOWe the c. T688211 MOW the C 133 aliquot of cells that had been growing in the shake flask for 22 hours was inoculated in the center of the chamber. Because of a problem with the camera, no images were taken until 16 hours after inoculation. At that time, a very unique pattern had formed in the DGC, as shown in Figure 57A. This pattern includes spots indicating local high concentrations of E. coli. The pattern continued to evolve, as shown in Figure 57B (18 hours) and Figure 57C (20 hours). Attempts to reproduce these results failed. Figure 57. Pattern formation in the DGC. A version of the mathematical was developed to incorporate a chemoattractant produced by the bacteria. This model was tested, and appeared to be working properly. However, no patterns were observed in the simulations. Work may be necessary to find the correct parameter combinations that allow for pattern formation, if this line of research is continued. Movement around objects A group of swarm plate experiments was performed to analyze the movement of the chemotactic wave of Pseudomonas KC around objects. The experimental conditions are des object channe around was [11: the cor cells vs gel sur hypotl apprm be har in ca; capillz 134 are described in Section 7.1.3. In Figure 58, three objects were placed in the gel. The object in the lower right—hand section of the gel was a tubing connector with an open channel through it. It appeared that the cells swam more quickly through the opening than around the object. Two hypotheses were formed to explain this phenomenon. The first was that the tubing connector was full of water remaining from the autoclave cycle. When the connector was immersed in the gel, the gel may not have displaced the water. The cells would then encounter less resistance to swimming in the water as compared to the gel surrounding the connector, and could move through the connector faster. The second hypothesis was that, in the narrow confines of the connector (the channel is approximately 0.4 mm in diameter), the cells tumbling would be confined, and it would be harder for the cells to change swimming direction. Evidence of this has been observed in capillary assays (Liu and Papadopoulos, 1995; Liu et al., 1996), although their capillaries had much smaller diameters. Figure 58. Chemotaxis through a hollow object. APPENDIX B APPEl The All and fro tndiagc time st: APPENDIX B The Alternating-Direction Implicit Method This discussion and the figures have been adapted from Chapra and Canale (1988) and from Carnahan et al. ( 1969). The altemating—direction implicit (ADD scheme uses tridiagonal matrices to solve parabolic equations in two or more spatial dimensions. Each time step is broken into two steps, as shown in . o Explicit o Implicit First half-step Second half-step Figure 59. The alternating-direction implicit ethod. The first step to solve the simple two-dimensional diffusion equation given in Equation ( 4) is written 135 such tl implici Equatiu where tridiag apprO) so tha explic Again 136 —1.j 1,j+l .,j—1 At / 2 ’ 5 (M2 + (Ay)2 SSW — Silij D 51:14 ‘25:]. + S," Sim/2 - 25:?” +S."+l/2] (53) such that the approximation of 828/ax2 is explicit, while the approximation of 828/8y2 is implicit. For the model presented in Section 4.2, the grid was square, and thus Ax=Ay, so Equation ( 53) may be rearranged as — 1W2 + 2(1+ ”5,33%”? — 154“” = 15:11.]. + 2(1— A)S,'jj + 215,1, , ( 54) uj-l l,j+l where X=DSAU(Ax)2. When Equation ( 54) is written for each point on the grid, a tridiagonal set of simultaneous equations results. For the second half step, Equation ( 4) is approximated by Stiff _ Si’fil/z _ Sir-ii; _ 25:31 + Sin—iij + Still/12 — 25:?” + 5:;32 ( 55) At / 2 5 (my (Ayf so that now the approximation of azs/axz is implicit, and the approximation of 828/8y2 is explicit. Equation ( 55) can be rearranged to yield — 15:17}, + 2(1+ MS“ — 15,17} , = 2.53”“ + 2(1— ”5,7?” + 25,737,112 ( 56) ivj [hi—1 Again, when Equation ( 56) is written for the entire grid, a tridiagonal matrix results. APPENDIX C APPE Instrur Th (ua ant chamo that th< 0 Th 0 Fil OEa Onaf diffus Spaci the g grow eXpe outp APPENDIX C Instructions for FORTRAN model The program currently solves five coupled PDE'S. These include two cell balances (ua and ub), two chemoattractants (s and q) and a growth nutrient (h). The sides of the chamber are identified by compass directions (n,s,e,w). For the version of the program that these instructions pertain to, the chemoattractant q should be thought of as oxygen. 0 The unit basis for all variables is mass=grams, length=centimeters, time=hours. 0 Filenames are italicized, while variable names are in "quotes". 0 Each input file is formatted, so be sure to keep numbers in the proper positions. On a Sun or HP: 1. Modify the param file. This file contains the parameters of the model, including chemotactic sensitivities, diffusion coefficients, mass transport coefficients, etc. The param file also contains information about the simulation run-time and grid Spacing. The value of "m" sets the grid size (m x m matrix). "Timeinit" sets the time that the gradients establish before cell inoculation. "Timeinoc" sets the time that the cells grow in the chamber. "Timeinit" plus "timeinoc" gives the total simulated time of the experiment. Param also contains instructions for the simulation that specify a movie output or single graph output. 137 _' chemt begin: exam] gel is begin: conce initial and "1 cell in "Num four g "time: Para" not pr 4. n S “36d. 138 A list of variables contained in the param file and their definition is included at the end of this section. 2. Modify the init file. This file sets the initial conditions, including initial cell concentration, chemoattractant boundary conditions, and nutrient boundary conditions. The parameters beginning with r (rsn, for example) correspond to the reservoir concentrations. In the example, rsn is reservoir, s-chemoattractant, north side. The initial concentration in the gel is also set in init. For example, "stini" sets the concentration of s in the gel at the beginning of the simulation. The nutrient concentration is usually set equal to the nutrient concentration in the four reservoirs. For a chemoattractant that will make a gradient, the initial concentration in the gel is usually 0. The cell concentration is initialized by "uaint" and "ubini". Note that these are not the initial concentrations, but are a parameter for the cell initialization subroutine. 3. Modify the "time" file. "Numti" tells the program how many matrices to print. For example, for "numti" =04, four graphs would be produced at times specified by "timel "time2", etc. The times are time after inoculation. Note that the variable "timeinoc" in param specifies how long the simulation will run, so any time longer than "timeinoc" will not print a graph. 4. If necessary, modify the FORTRAN program itself. Subroutine uaO sets the initial cell distribution. Currently, an exponential peak is used. Tl both f be cha differe Ti 0 bound comp‘ Toru word Dati The Calle Whh 139 The subroutines f and ifs set the nutrient and chemoattractant uptake functions. Note: both f and fs must have identical functionalities. If one is changed, the other must also be changed. Both functions correspond to the same term in the equations, but are slightly different because of the finite differencing algorithm. The subroutine chemofl sets the functionality of the chemotaxis term. Other subroutines set the windows of the reservoirs, calculate velocities, set boundary conditions, and calculate certain other values. 5. To run the program, a Sun computer must be used. If an HP is the current computer, open a window to a Sun. To run the program in batch mode, type: batch An at> prompt should appear. Type the following commands (..J = return, underlined words should be replaced with appropriate name, and bold indicates pressing control and D at the same time): at> f77 _tllegamej .J at> a.out .J at> mailx -s "subject" mail-address J at> ctrl-D The simulation is complete when mail bearing the "subject" is received. At this time, files called lcella.m, lcellb.m, latts.m, lattq.m, and lnut.m should be in the directory from which the program ran. These files are formatted for Matlab. Matlal 1.A 2.A 3. Ty exam 4. I For 6 shoul 5. Tv will Higl give Mat help 140 Matlab instructions: 1. At a Sun or HP, start Matlab by typing: matlab 2. At the Matlab prompt, switch to the directory containing the simulation files: cd path\directog 3. Type the name of the file you wish to work with, without the m extension. For example, to look at the cell Population a matrix, type lcella 4. To check if the matrices have been properly loaded, type who For example, if "numti" = 4, and if individual graphs were specified, then there should be 4 matrices named ual, ua2, ua3, and ua4. 5. Two types of graphs are useful for viewing: Three dimensional plots are created with the surf command. surf(ual) will yield a 3-D plot, with Matlab default view angle and axis. Overhead view type plots are created with the appmult command. appmult(ual,3,'pcolor') will give a good overhead plot. The value 3 can be changed to give better or worse plots. Higher numbers give smoother graphs, but take longer to view and print. Lower numbers give coarser plots, but are faster. Matlab has a fairly good help feature. Typing help will list general help topics, or type help command to get specific help. APPENDIX D APP: Thel canb Popu APPENDIX D The FORTRAN model This is the FORTRAN model used to solve the competition model. This model can be reduced to the one-cell population model simply by setting the initial condition for Population B to zero in the initialization file (see Appendix E). implicit double precision (a—h,o-z) C parameter (nmax=100,kmax=100) C dimension x(0:nmax),y(0:nmax) dimension ual(O:nmax,0:nmax),ua2(0:nmax,0:nmax) dimension ub1(0:nmax,0:nmax),ub2(0:nmax,0:nmax) dimension H1(O:nmax,0:nmax),H2(O:nmax,O:nmax) dimension S 1(O:nmax,0:nmax),82(0:nmax,0:nmax) dimension q1(O:nmax,0:nmax),q2(0:nmax,0:nmax) dimension aa(0:nmax),bb(0:nmax),cc(0:nmax) dimension funua(O:nmax),funs(0:nmax),funh(0:nmax) dimension rhsua(0:nmax),rhss(0:nmax),rhsh(0:nmax) dimension funub(0:nmax) dimension rhsq(0:nmax),funq(0:nmax),rhsub(0:nmax) dimension uua(0:nmax),us(0:nmax),uh(0:nmax),uq(0:nmax),uub(0:nmax) dimension mcount3(0:kmax), time(Ozkmax) dimension velm(0:500),vell(0:500) dimension distm(0:500),distl(0:500),timm(0:500),timl(0:500) open(16,file='lcella.m') open(32,file='lcellb.m’) open( 17,file='latts.m') open(23 ,file='1attq.m') open(18,file='lnut.m') open( 1 9,fi1e='time') open(24,file=‘init') open(22,file=‘param') open(30,file='distmax') open(3 1 ,file='distlim') 141 142 read(22,101)m read(22, 103)timinit read(22,103)timinoc read(22, 102)numcel read(22,102)numvel read(22,102)kstep read(22,100)dmin read(22,102)movc read(22,100)R read(22, 100)ayoff read(22, 100)axoff read(22, 100)byoff read(22, 100)bxoff read(22, 100)awid read(22, 100)bwid read(22,100)Dh read(22, 100)thn read(22, 100)ths read(22, 100)the read(22, 100)thw read(22, 100)Ds read(22, 100)tsn read(22, 100)tss read(22, 100)tse read(22, 100)tsw read(22,100)Dacs read(22, 100)Dbcs read(22, 100)DKas read(22, 100)DKbs read(22, 100)cas read(22, 100)vas read(22, 100)Yas read(22, 100)cbs read(22,100)vbs read(22,100)Ybs read(22, 100)Dq read(22, 100)tqn read(22, 100)th read(22, 100)tqe read(22, 100)tqw read(22,100)Dacq read(22, 100)Dbcq read(22, 100)DKaq 143 read(22, 100)DKbq read(22, 100)caq read(22, 100)vaq read(22, 100)Yaq read(22, 100)cbq read(22, 100)qu read(22, 100)qu read(22,100)Dua read(22,100)va read(22,100)ca read(22, 100)Ya read(22,100)Dub read(22, 100)vb read(22, 100)cb read(22,100)Yb read(24, 100)uaini read(24,100)ubini read(24, 100)rsn read(24,100)rss read(24, 100)rse read(24,100)rsw read(24, 100)rqn read(24, 100)th read(24, 100)rqe read(24,100)rqw read(24, 100)rhn read(24, 100)rhs read(24, 100)rhe read(24, 100)rhw read(24, 100)stini read(24, 100)qtini read(24, 100)htini 100 format(6x,dl3.10) 101 format(4x,i4) 102 format(7x,i7) 103 format(8x,d13.10) if(movc.eq.2)then write(l6,*)'ua=[' write(32,*)'ub=[' write(17,*)'s=[' write(23,*)'q=[' write(18,*)'h=[' endif C 144 hr=2.0d0*R/m ht=0.15625d0*hr*hr mcount1=timinit/ht+1 mcount2=timinocfht+1 write(*,*)'hr =',hr write(*,*)'ht =',ht write(*,*)'mcount1=',mcount1 write(*,*)'mcount2=',mcount2 read(19,201)numti do 555 i=1,numti read(19,202)time(i) mcount3(i)=time(i)/ht time(i)=mcount3(i)*ht write(*,*)'Real time',i,'=',time(i) 555 continue 201 forrnat(6x,i2) 202 format(6x,d6.2) C Calculate values of hr and ht C C set up the initial gradient of attractant and nutrient C C set up the initial chemoattractant concentration in chamber do 10 i=0,m do 5 j=0,m S l (i,j)=stini q1(i,j)=qtini H1(i,j)=htini ua1(i,j)=0.0d0 ub1(i,j)=0.0d0 5 continue 10 continue C D0=Dh*(ht/2.0d0)/(hr**2) D1=Ds*(ht/2.0d0)/(hr**2) D2=Dua*(ht/2.0d0)/(hr**2) D3=Dacs*(ht/2.0d0)/(hr**2) D4=Dq*(ht/2.0d0)/(hr**2) D5=Dacq*(ht/2.0d0)/(hr**2) D6=Dub*(ht/2.0d0)/(hr**2) D7=Dbcs*(ht/2.0d0)/(hr**2) D8=Dbcq*(ht/2.0d0)/(hr**2) C C solving the ode 145 do 1000 k=1,mcount1 C C take care of boundary points: C i=0 do 20 i=0,m t=tr(i,hr,tss,Ds) call rhbound(m,S1(i,0),S1(i,1),rss,D1,hr,t,rhss(i)) funs(i)=1.0d0 t=tr(i,hr,tqs,Dq) call rhbound(m,q1(i,0),q1(i,1),rqs,D4,hr,t,rhsq(i)) funq(i)=1 .0d0 t=tr(i,hr,ths,Dh) call rhbound(m,H1(i,O),H1(i,1),rhs,D0,hr,t,rhsh(i)) funh(i)=1 .0d0 20 continue C call PDMTRIX(m,D1,funs,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhss,us) call PDMTRIX(m,D4,funq,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhsq,uq) call PDMTRIX(m,DO,funh,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhsh,uh) do 25 i=0,m 82(i,0)=us(i) q2(i.0>=uq=82(i.j> q1(1.1)=q2=uho> ua2(0,j)=uua(j) ub2(0,j)=uub(j) continue i=m rhss(0)=Sl(m,0)+D1*(S 1(m—1,0)-S 1(m,0)) rhss(0)=rhss(0)-(ht/2.0d0)*f(cas,vas,s 1(m,0))*ua1 (m,0)/Y as rhss(0)=rhss(0)-(ht/2.0d0)*f(cbs,vbs,s1(m,0))*ub1(m,O)/Ybs funs(0)=1 .0d0 rhsq(0)=q1(m,0)+D4*(q1(m-1,0)—q1(m,0)) rhsq(0)=rhsq(0)-(ht/2.0d0)*f(caq,vaq,q1(m,0))*ua1(m,0)/Yaq rhsq(0)=rhsq(0)—(ht/2.0d0)*f(cbq,qu,q1(m,0))*ub1(m,0)leq funq(0)=1 .OdO rhsh(0)=h1(m,0)+DO*(h1(m-1,0)-h1(m,0)) rhsh(0)=rhsh(0)—(ht/2.0d0)*f(ca,va,h1(m,0))*ua1(m,0)/Ya rhsh(0)=rhsh(0)-(ht/2.0d0)*f(cb,vb,h1(m,0))*ub1(m,0)/Yb funh(0)=1 .0d0 rhss(m)=S 1(m,m)+D1 *(S 1(m-1,m)—S 1(m,m)) rhss(m)=rhss(m)~(ht/2.0d0)*f(cas,vas,s 1 (m,m))*ual (m,m)/Y as rhss(m)=rhss(m)—(ht/2.0d0)*f(cbs,vbs,s1(m,m))*ub1(m,m)/Ybs funs(m)=1.0d0 rhsq(m)=q1(m,m)+D4*(ql(m-1,m)-q1(m,m)) rhsq(m)=rhsq(m)-(ht/2.0d0)*f(caq,vaq,q1(m,m))*ua1(m,m)/Yaq rhsq(m)=rhsq(m)-(ht/2.0d0)*f(cbq,qu,q1(m,m))*ub1(m,m)/qu funq(m)=1 .0d0 rhsh(m)=hl(m,m)+D0*(h1(m-1,m)—h1(m,m)) rhsh(m)=rhsh(m)—(ht/2.0d0)*f(ca,va,h1(m,m))*ual (m,m)/Ya rhsh(m)=rhsh(m)-(ht/2.0d0)*f(cb,vb,hl (m,m))*ub l(m,m)/Yb funh(m)=1 .0d0 C Cell group a a21=0.5d0*chemof1(DKas,S1(m-1,0))*ua1(m-1,0) a21=a2l+0.5d0*chemof1(DKas,S 1(m,0))*ual(m,0) b12=0.5d0*chemof1(DKas,S1(m,1))*ua1(m,1) b12=bl2+0.5d0*chemof1(DKas,S 1(m,0))*ua1(m,0) 021:0.5d0*chemof1(DKaq,q1(m—1,0))*ua1(m—1,0) c21:021+0.5d0*chemof1(DKaq,q1(m,0))*ua1(m,0) d12=0.5d0*chemof1(DKaq,q1 (m, 1))*ua1(m, 1) d12=d12+0.5d0*chemof1(DKaq,q1(m,0))*ua1(m,0) rhsua(O)=ua1(m,0)+D2*(ua1(m-1,0)-ua1(m,0)) rhsua(O)=rhsua(0)+(ht/2.0d0)*f(ca,va,h1(m,0))*ua1(m,0) rhsua(O)=rhsua(0)—D3 *a21 *(s 1(m— 1 ,0)-sl(m,0)) rhsua(O)=rhsua(0)-D3 *b12*(S 1(m, 1)-s 1(m,0)) 170 rhsua(O)=rhsua(0)-D5 *c21 *(q 1 (m- 1,0)-q1(m,0)) rhsua(O)=rhsua(0)-D5 *d 1 2* (q 1 (m, 1)-q1 (m,0)) funua(0)= 1 .OdO a21=0.5d0*chemof1(DKas,S1(m-1,m))*ua1(m—1,m) a21=a21+0.5d0*chemof1(DKas,S1(m,m))*ua1(m,m) b21=0.5d0*chemofl (DKas,S 1(m,m— l ))*ual (m,m— 1) b21=b21+0.5d0*chemof1(DKas,S 1(m,m))*ua1(m,m) c2l=0.5d0*chemof1(DKaq,ql(m-1,m))*ual(m-1,m) c21=c21+0.5d0*chemof1(DKaq,q1(m,m))*ual(m,m) d21=0.5d0*chemof1(DKaq,q1(m,m-1))*ua1(m,m-1) d21=d21+0.5d0*chemof1(DKaq,q1(m,m))*ua1(m,m) rhsua(m)=ua1(m,m)+D2*(ua1(m-1,m)—ua1(m,m)) rhsua(m)=rhsua(m)+(ht/2.0d0)*f(ca,va,h1(m,m))*ua1(m,m) rhsua(m)=rhsua(m)-D3 *a21 *(sl(m-1,m)-sl(m,m)) rhsua(m)=rhsua(m)-D3 *b21*(sl(m,m—1)—sl(m,m)) rhsua(m)=rhsua(m)-D5 *c21*(q1(m—1,m)-q1(m,m)) rhsua(m)=rhsua(m)—D5 *d21*(q1(m,m—l)—q1(m,m)) funua(m)=1.0d0 C C Cell group b a21=0.5d0*chemof1(DKbs,S 1(m-1,0))*ub1(m—1,0) a21=a21+0.5d0*chemof1(DKbs,S 1(m,0))*ub1(m,0) b 12=0.5d0*chemof1(DKbs,S 1(m, 1))*ub1(m, 1) b12=b12+0.5d0*chemof1(DKbs,S1(m,0))*ub1(m,0) 021=0.5d0*chemof1(DKbq,ql(m-1,0))*ub1(m—1,0) 021:021+0.5d0*chemof1(DKbq,q1(m,0))*ub1(m,0) d12=0.5d0*chemof1(DKbq,q1(m, 1))*ub1(m,1) d12=d12+0.5d0*chemof1(DKbq,ql(m,0))*ub1(m,0) rhsub(0)=ub1(m,0)+D6*(ub](m—1,0)-ub1(m,0)) rhsub(0)=rhsub(0)+(ht/2.0d0)*f(cb,vb,h1(m,0))*ub1(m,0) rhsub(0)=rhsub(0)-D7*a21 *(s 1(m—l ,0)-s 1(m,0)) rhsub(0)=rhsub(0)-D7*b12*(S 1(m, 1)-s 1(m,0)) rhsub(0)=rhsub(0)-D8*c21*(ql(m—1,0)-q1(m,0)) rhsub(0)=rhsub(0)-D8 *d12*(q1(m,1)-q1(m,0)) funub(0)=-1.0d0 a21=0.5d0*chemof1(DKbs,S1(m-1,m))*ub1(m-1,m) a21=a21+0.5d0*chemof1(DKbs,S1(m,m))*ubl(m,m) b21 =0.5d0*chemof1(DKbs,S 1(m,m— 1))*ub1(m,m- 1) b21=b21+0.5d0*chemof1(DKbs,S l(m,m))*ub1(m,m) c21=0.5d0*chemof1(DKbq,q1(m-1,m))*ubl(m—1,m) 171 c21=c2l+0.5d0*chemof1(DKbq,q1(m,m))*ub1(m,m) d21=0.5d0*chemof1(DKbq,q1(m,m-1))*ub1(m,m—1) d21=d21+0.5d0*chemof1(DKbq,ql(m,m))*ub1(m,m) rhsub(m)=ub1(m,m)+D6*(ub1(m-1 ,m)-ub1(m,m)) rhsub(m)=rhsub(m)+(ht/2.0d0)*f(cb,vb,h1(m,m))*ub 1(m,m) rhsub(m)=rhsub(m)—D7*a21*(s1(m-1,m)-s1(m,m)) rhsub(m)=rhsub(m)-D7*b21 *(sl(m,m—1)-sl(m,m)) rhsub(m)=rhsub(m)-D8*c21*(ql(m—1,m)-q1(m,m)) rhsub(m)=rhsub(m)-D8 *d21*(q1(m,m- 1)-q 1 (m,m)) funub(m)=1.0d0 do 1070 j=1,m-1 t=tr(i ,hr,tse,Ds) call rhbound(m,S 1 (m,j),S 1(m-1,j),rse,D1,hr,t,rhss(i)) rhss(i)=rhss(j)-(ht/2.0d0)*f(cas,vas,s1(m,j))*ua1(m,j)/Y as rhss(i)=rhss(i)-(ht/2.0d0)*f(cbs,vbs,s1(m,j))*ub1(m,j)/Ybs funs(0)=1 .0d0 t=tr(j ,hr,tqe,Dq) call rhbound(m,ql(m,j),ql(m—1,j),rqe,D4,hr,t,rhsq(i)) rhsq(j)=rhsqG)-(ht/2.0d0)*f(caq,vaq,q1(m,j))*ua1(m,j)/Yaq rhsq(j)=rhsq(j)-(ht/2.0d0)*f(cbq,qu,q1(m,j))*ub1(m,j)/qu funq(0)=1 .0d0 t=tr(j ,hr,the,Dh) call rhbound(m,hl(m,j),h1(m—1,j),rhe,D0,hr,t,rhsh(j)) rhsh(j)=rhsh(j)-(ht/2.0d0)*f(ca,va,h1(m,j))*ua1(m,j)/Ya rhsh(j)=rhsh(j)-(ht/2.0d0)*f(cb,vb,h1(m,j))*ub1(m,j)/Yb funh(0)=1 .0d0 C Cell group a a21=0.5d0*chemof1(DKas,S1(m-1,j))*ua1(m—1,j) a21=a21+0.5d0*chemof1(DKas,S1(m,j))*ua1(m,j) b12=0.5d0*chemof1(DKas,S1(m,j+1))*ual(m,j+1) b12=b12+0.5d0*chemof1(DKas,S1(m,j))*ua1(m,j) b21=0.5d0*chemof 1(DKas,S 1(m,j- 1))*ua1(m,j- 1) b21=b21+O.5d0*chemof1(DKas,S1(m,j))*ua1(m,j) c21=0.5d0*chemof1(DKaq,q1(m-1,j))*ua1(m-1,1) c21=c2l+0.5d0*chemof1(DKaq,ql (m,j))*ua1(m,j) d12:0.5d0*chemof1(DKaq,ql(m,j+1))*ua1(m,j+1) d12=d12+0.5d0*chemof1(DKaq,ql (m,j))*ua1(m,j) d21=0.5d0*chemof1(DKaq,q1(m,j-1))*ua1(m,j-1) d21=d21+0.5d0*chemof1(DKaq,ql(m,j))*ua1(m,j) rhsua(i)=ua](m,j)+D2*(ua1(m-1,j)—ua1(m,j)) rhsua(i)=rhsua(j )+(ht/2.0d0)*f(ca,va,h1(m,j))*ua1(m,j) rhsua(i)=rhsua(j)-D3*a21 *(sl(m—1,j)-sl(m,j)) 172 rhsua(i)=rhsua(j)-D3*b12*(S1(m,j+1)-s1(m,j)) rhsuafi)=rhsua(j)-D3*b21*(S1(m,j-1)-sl(m,j)) rhsua(i)=rhsuaG)-D5*c21*(ql(m—1,j)-q1(m,j)) rhsua(i)=rhsua(i)-D5*d12*(q1(m,j+1)-q1(m,j)) rhsua(i)=rhsua(j)-D5*d21*(ql(m,j-1)-ql(m,j)) funua(i)=1 .0d0 C C Cell group b a21=0.5d0*chemof1(DKbs,Sl(m—1,j))*ub1(m-1,j) a21=a21+0.5d0*chemof1(DKbs,S1(m,j))*ub1(m,j) b12=0.5d0*chemof1(DKbs,S l(m,j+1))*ub1(m,j+1) b12=b12+0.5d0*chemof1(DKbs,S1(m,j))*ub1(m,j) b21=0.5d0*chemofl(DKbs,S 1(m,j-1))*ub 1(m,j-1) b21=b21+0.5d0*chemof1(DKbs,S1(m,j))*ub1(m,j) c21=0.5d0*chemof1(DKbq,ql(m-1 ,j))*ub1(m-1,j) c21=c21+0.5d0*chemof1(DKbq,q1(m,j))*ubl (m,j) d12=0.5d0*chemof1(DKbq,q1(m,j+1))*ub1(m,l+1) d12=d12+0.5d0*chemof1(DKbq,q1 (m,j))*ub 1(m,j) d21=0.5d0*chemof1(DKbq,q1(m,j-1))*ub1(m,j-1) d21=d2 1+0.5d0*chemof1 (DKbq,q1(m,j ))*ub1(m,j) rhsub(i)=ub1(m,j)+D6*(ub1(m-1,j)—ub1(m,j)) rhsub(i)=rhsub(j )+(ht/2.0d0)*f(cb,vb,h1 (m,j)) *ub 1 (m, j) rhsub(j )=rhsub(j)-D7 *a2 1 *(s 1 (m— 1 ,j)-s 1 (m,j )) rhsub(j)=rhsub(j)-D7*b12*(S 1(m,j+1)-s 1 (m,j)) rhsub(i)=rhsub(i)-D7*b21*(S1(m,j-1)-sl(m,j)) rhsub(j)=rhsub(j)-D8*c21*(ql(m-l,j)-q1(m,j)) rhsub(i)=rhsub(i)—D8*d12*(q1(m,j+1)—q1(m,j)) rhsub(i)=rhsub(j)-D8*d21*(ql(m,j-l)-q1(m,j)) funub(i)=1 .0d0 1070 continue C call PDMTRIX(m,D1,funs,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhss,us) call PDMTRIX(m,D4,funq,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhsq,uq) call PDMTRD((m,DO,funh,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhsh,uh) ' call PDMTRlX(m,D2,funua,aa,bb,cc) call TRIS OLV(m,aa,bb,cc,rhsua,uua) call PDMTRIX(m,D6,funub,aa,bb,cc) call TRISOLV(m,aa‘,bb,cc,rhsub,uub) do 1075 j=0,m sztm.j)=us0) q2(m,j)=uq(i) 173 h2(m,j)=uh(i) ua2(m,j)=uua(j) ub2(m,j)=uub0') 1075 continue C take care of interior points: C do 1090 i=1,m- 1 C t=tr(i,hr,tss,Ds) call rh2(m,S1(i,0),Sl(i—1,0),S1(i+1,0),rss,D1,hr,t,rhss(0)) rhss(O)=rhss(0)-(ht/2.0d0)*f(cas,vas,s1(i,0))*ua1(i,0)/Yas rhss(0)=rhss(0)-(ht/2.0d0)*f(cbs,vbs,s1(i,0))*ub1 (i,0)/Ybs funs(0)=l .0d0+D1*hr*tr(i,hr,tss,Ds) t=tr(i,hr,tqs,Dq) call rh2(m,q1(i,0),q1(i-1,0),q1(i+1,0),rqs,D4,hr,t,rhsq(0)) rhsq(0)=rhsq(O)-(ht/2.0d0)*f(caq,vaq,q1(i,0))*ua1(i,0)/Yaq rhsq(0)=rhsq(0)-(ht/2.0d0)*f(cbq,qu,q1(i,0))*ub1(i,0)/qu funq(0)=1.0d0+D4*hr*tr(i,hr,tqs,Dq) t=tr(i,hr,ths,Dh) call rh2(m,h1(i,0),h1(i-1,0),h1(i+1,0),rhs,D0,hr,t,rhsh(0)) rhsh(0)=rhsh(0)-(ht/2.0d0)*f(ca,va,h1(i,0))*ua1(i,0)/Ya rhsh(0)=rhsh(0)-(ht/2.0dO)*f(cb,vb,h1(i,0))*ub1(i,0)/Yb funh(0)=1 .0d0+D0*hr*tr(i,hr,ths,Dh) t=tr(i,hr,tsn,Ds) call rh2(m,S1(i,m),S1(i-1,m),S1(i+1,m),rsn,D1,hr,t,rhss(m)) rhss(m)=rhss(m)-(ht/2.0dO)*f(cas,vas,s1(i,m))*ua1(i,m)/Yas rhss(m)=rhss(m)-(ht/2.0d0)*f(cbs,vbs,s1(i,m))*ub1(i,m)/Ybs funs(m)=1.0d0+D1*hr*tr(i,hr,tsn,Ds) t=tr(i,hr,tqn,Dq) call rh2(m,ql(i,m),q1(i-1,m),q1(i+1,m),rqn,D4,hr,t,rhsq(m)) rhsq(m)=rhsq(m)—(ht/2.0d0)*f(caq,vaq,q1(i,m))*ua1(i,m)/Yaq rhsq(m)=rhsq(m)-(ht/2.0d0)*f(cbq,qu,q1(i,m))*ub1(i,m)/qu funq(m)=1 .0d0+D4*hr*tr(i,hr,tqn,Dq) t=tr(i,hr,thn,Dh) call rh2(m,h1(i,m),h1(i-1,m),hl(i+1,m),rhn,D0,hr,t,rhsh(m)) rhsh(m)=rhsh(m)—(ht/2.0dO)*f(ca,va,hl(i,m))*ua1(i,m)/Ya rhsh(m)=rhsh(m)—(ht/2.0d0)*f(cb,vb,h1 (i,m))*ub1(i,m)/Yb funh(m)=1 .0d0+D0*hr*tr(i,hr,thn,Dh) C i=0 C Cell group a a1 2=0.5d0*chemof1(DKas,S 1(i+1 ,0))*ua1(i+1,0) a12=a12+0.5d0*chemof1(DKas,S1(i,0))*ua1(i,0) ('30 174 a21=0.5d0*chemof1(DKas,S1(i,0))*ua1(i,0) a21=a21+0.5d0*chemof1(DKas,S 1(i-1,0))*ua1(i- 1,0) b12=0.5d0*chemof1(DKas,S 1(i,1))*ua1 (i, l) b12=b12+0.5d0*chemof1(DKas,S1(i,0))*ual(i,0) c12=0.5d0*chemof1(DKaq,q1(i+1,0))*ua1(i+1,0) c 12=c12+0.5d0*chemof1(DKaq,q1(i,0))*ua1(i,0) c21=0.5d0*chemof1(DKaq,q1(i,0))*ua1(i,0) 021=c21+0.5d0*chemof1(DKaq,q1(i-1 ,0))*ua1(i- 1 ,0) d12=0.5d0*chemof1(DKaq,q](i,1))*ua1(i,1) d12=d12+0.5d0*chemof1(DKaq,q1(i,0))*ua1(i,0) . rhsua(O)=ua1(i,0)+D2*(ua1(i+1,0)-ua1(i,0)) * +D2*(ual (i-1,0)-ua1(i,0)) rhsua(O)=rhsua(0)+(ht/2.0d0)*f(ca,va,h1(i,0))*ua1(i,0) rhsua(O)=rhsua(0)-D3 *a12*(s 1(i+1,0)-s 1(i,0)) rhsua(O)=rhsua(0)-D3*a21*(sl(i-1,0)-sl(i,0)) ' rhsua(O)=rhsua(0)-D3*b12*(S1(i,1)-sl(i,0)) rhsua(O)=rhsua(0)-D5*c12*(q1(i+1,0)-q1(i,0)) rhsua(O)=rhsua(O)-D5*021*(q1(i-1,0)-q1(i,0)) rhsua(O)=rhsua(0)—D5 *d12*(q1(i,1)-q1(i,0)) funua(0)=l .0d0 i=m a12=0.5d0*chemofl(DKas,S l(i+1,m))*ua1(i+1,m) a12=a12+0.5d0*chemof1(DKas,S1(i,m))*ua1(i,m) a21=0.5d0*chemof1(DKas,S1(i,m))*ua1(i,m) a21=a21+0.5d0*chemof1(DKas,S 1(i-1,m))*ua1(i-1,m) b21=0.5d0*chemofl(DKas,S1(i,m))*ua1(i,m) b21=b21+0.5d0*chemof1(DKas,S1(i,m-1))*ua1(i,m-1) c12=0.5d0*chemof1(DKaq,q1(i+1,m))*ua1(i+1,m) c12=c12+0.5d0*chemof1(DKaq,q1(i,m))*ua1(i,m) c21=0.5d0*chemof1(DKaq,ql(i,m))*ual(i,m) c2l=c21+0.5d0*chemof1(DKaq,ql(i-1 ,m))*ual(i-1 ,m) d21=0.5d0*chemof1(DKaq,ql (i,m))*ua1(i,m) d21=d2l+0.5d0*chemof1(DKaq,q1(i,m-1))*ua1(i,m-1) rhsua(m)=ua1(i,m)+D2*(ua1(i+1,m)-ual(i,m)) * +D2*(ua1(i-1,m)-ual(i,m)) rhsua(m)=rhsua(m)+(ht/2.0d0)*f(ca,va,h1(i,m))*ua1(i,m) rhsua(m)=rhsua(m)—D3 *al2*(s1(i+1,m)-Sl(1,m)) rhsua(m)=rhsua(m)-D3*a21*(sl(i-1,m)-Sl(i,m)) rhsua(m)=rhsua(m)-D3*b21*(S 1 (i,m-1)-s 1(i,m)) rhsua(m)=rhsua(m)-D5*c12*(q1(i+1,m)-q1(i,m)) rhsua(m)=rhsua(m)-D5*021*(ql(i-1,m)-q1(i,m)) C 175 rhsua(m)=rhsua(m)—D5*d21*(ql(i,m—1)-ql(i,m)) funua(m)= l .0d0 C Cell group b 00 a12=0.5d0*chemofl(DKbs,S1(i+1,0))*ub1(i+l,0) a12=312+0.5d0*chemof1(DKbs,Sl(i,0))*ub1(i,0) a21=0.5d0*chemof1(DKbs,S1(i,0))*ub1(1,0) a21=a21+0.5d0*chemof1(DKbs,S 1(1- 1 ,0))*ub1(i-1,0) b12=0.5d0*chemof1(DKbs,S1(i,1))*ub1(1, 1) b12=b12+0.5d0*chemof1(DKbs,S1(1,0))*ub1(i,0) c12=0.5d0*chemof1(DKbq,ql(i+1,0))*ub1(i+1,0) c12=c12+0.5d0*chemof1(DKbq,q1(i,0))*ubl(i,0) c21=0.5d0*chemof1(DKbq,q1(1,0))*ub1 (1,0) c21=c21+0.5d0*chemof1(DKbq,q1(i-1,0))*ub1(i—1,0) d12=0.5d0*chemof1(DKbq,ql(i, 1))*ub1(i, 1) d12=d12+0.5d0*chemof1(DKbq,q1(i,0))*ub 1 (1,0) rhsub(0)=ub1 (i,0)+D6*(ub 1(i+ 1 ,0)-uh 1(i,0)) +D6*(ub1(1—1,0)-ub1(i,0)) rhsub(0)=rhsub(0)+(ht/2.0d0)*f(cb,vb,h1 (i,0))*ub1 (1,0) rhsub(0)=rhsub(0)—D7*a12*(s1(i+1,0)-sl(i,0)) rhsub(0)=rhsub(0)—D7*a21 *(sl(i-1,0)—sl (1,0)) rhsub(0)=rhsub(0)—D7*b12*(S l (i,1)-sl(i,0)) rhsub(0)=rhsub(0)-D8*c12*(q1(i+1,0)—q1(i,0)) rhsub(0)=rhsub(0)—D8*c21*(ql(i-1,0)—q1(i,0)) rhsub(0)=rhsub(0)-D8*d12*(q1(i,1)—q1(i,0)) funub(0)=1.0d0 J=m a12=0.5d0*chemof 1 (DKbs,S1(i+1,m))*ubl(i+1,m) a12=a12+0.5d0*chemof1(DKbs,S1(i,m))*ub1(i,m) a21=0.5d0*chemof1(DKbs,Sl(i,m))*ub1(i,m) a21=a21+0.5d0*chem0f1(DKbs,S1(1-1,m))*ub1(i-1,m) b21=0.5d0*chemof1(DKbs,Sl(i,m))*ub1(i,m) b21=b2l+0.5d0*chem0f1(DKbs,S1(i,m-1))*ubl(i,m—1) c12=0.5d0*chemof1(DKbq,q1(i+1,m))*ub1(i+1,m) cl2=c12+0.5d0*chemof1(DKbq,q1(i,m))*ub1(i,rn) c21=0.5d0*chemof1(DKbq,q1(i,m))*ubl(i,m) c21=c21+0.5d0*chem0f1(DKbq,q1(i-1,m))*ub1(i—1,m) d21=0.5d0*chem0f1(DKbq,q1(i,m))*ubl(i,m) d21=d21+0.5d0*chern0f1 (DKbq,q 1 (i,m- 1 ))*ub l (i,m— 1) rhsub(m)=ub1(i,m)+D6*(ub1(i+1,In)-ub1(i,rn)) +D6*(ub1(i-1,m)-ubl(i,m)) 176 rhsub(m)=rhsub(m)+(ht/2.0d0)*f(cb,vb,h1(i,m))*ub 1 (i,m) rhsub(m)=rhsub(m)-D7*a12*(s1(1+1 ,m)—s 1 (i,m)) rhsub(m)=rhsub(m)—D7* a21 *(s 1 (1-1 ,m)-s1(i.m)) rhsub(m)=rhsub(m)-D7 *b21 *(S 1(1,m- l)-S 1(1,m)) rhsub(m)=rhsub(m)—D8*c12*(q1(i+1,m)-q1(i,m)) rhsub(m)=rhsub(m)—D8*c2 1 * (q 1 (1— 1 ,m)-q 1(1,m)) rhsub(m)=rhsub(m)-D8 *d21 *(q 1 (i,m- 1)-q1(i,m)) funub(m)= 1 .0d0 do 1080 j=1,m-1 rhss(j)=Sl(1,j)+D1*(S1(1+l,j)-2.0d0*Sl(i,j)+S1(i-1,j)) rhss(j)=rhss(i)-(ht/2.0d0)*f(cas,vas,s 1 (i,j))*ual (1, j)/Y as rhss(j)=rhss(j)-(ht/2.0d0)*f(cbs,vbs,s1(i,j))*ub1(i,j)/Ybs funs(i)=1.0d0 rhsq(j)=q1(1,j)+D4*(q1(1+1 ,j)-2.0d0*q1 (i,j)+q 1(1-1 ,j)) rhsq(j)=rhsq(j)-(ht/2.0d0)*f(caq,vaq,q1(i,j))*ual(i,j)/Yaq rhsq(j)=rhsq(j)—(ht/2.0d0)*f(cbq,qu,q1(i,j))*ub1(i,j)/qu funq(i)=1.0d0 rhsh(i)=h1(1,j)+D0*(h1(1+1,j)-2.0d0*h1(i,j)+h1(i-1,j)) rhsh(i)=rhsh(j)-(ht/2.0d0)*f(ca,va,h1(i,j))*ual (i,j )/Y a rhsh(i)=rhsh(j)-(ht/2.0d0)*f(cb,vb,h1(i,j))*ub1(i,j)/Yb funh(i)= 1 .OdO C Cell group a al2=0.5d0*chemof1(DKas,S1(1+1,j))*ua1(i+1,j) a12=a12+0.5d0*chemof1(DKas,S1(i,j))*ua1(1,j) a21=0.5d0*chemof1(DKas,S 1(i,j))*ual(1,j) a21=a21+0.5d0*chemof1(DKas,S1(1-1,j))*ua1(1—1,j) b12=0.5d0*chemof1(DKas,S l (1,j+1))*ua1(1,j+1) b12=b12+0.5d0*chemof1(DKas,S1(1,j))*ua1(1,j) b21=0.5d0*chemofl(DKas,S1(i,j))*ua1(i,j) b21=b21+0.5d0*chemof1(DKas,S 1(1,j- 1))*ua1(i,j-1) cl2=0.5d0*chemofl(DKaq,ql(1+1,j))*ua1(i+1,j) c 12=c12+0.5d0*chemof1(DKaq,ql(i,j))*ua1 (i,j) c21=0.5d0*chemof1(DKaq,ql(i,j))*ua1(i,j) c21=c21+0.5d0*chemof1(DKaq,ql(i-1,j))*ua1(i-1,j) d12=0.5d0*chemof1(DKaq,q1(i,j+1))*ua1(i,j+1) d12=d12+0.5d0*chemof1(DKaq,ql(1,j))*ua1(1,j) d21=0.5d0*chemof1(DKaq,ql(i,j))*ual(i,j) d21=d21+0.5d0*chemof1(DKaq,ql(1,j-1))*ua1(i,j—1) rhsua(i)=ua1(i,j)+D2*(ua1(i+1,j)—ua1(i,j)) * +D2*(ua1(i-1,j)-ua1(1,j)) rhsua(i)=rhsua(j)+(ht/2.0d0)*f(ca,va,h1(i,j))*ua1(i,j) rhsua(i)=rhsua(j)-D3*a12*(sl(1+1,j)-s1(i,j)) 177 rhsua(i)=rhsua(i)-D3*a21*(s1(1- 1 ,j)—sl(i,j)) rhsua(i)=rhsua(i)-D3*b12*(S1(i,j+1)-S1(i,j» rhsua(i)=rhsua(i)—D3 *b2 1 *(s 1 (i,j-1)-S 1 (i,j)) rhsua(i)=rhsua(j)-D5*c12*(q1(i+1,j)-q1(1,j)) rhsan')=rhsua(i)—D5*021*(q1 (1—1 ,j)—ql (i,j)) rhsua(i)=rhsua(j)—D5*d12*(q1(i,j+1)-(110,11) rhsua(i)=rhsua(j)—D5*d21*(ql(1,j-l)-q1(i,j)) funua(i)=1 .OdO C Cell group b a12=0.5d0*chemofl(DKbs,S1(1+1,j))*ub1(i+1,j) a12=a12+0.5d0*chemof1(DKbs,S1(i,j))*ub1(i,j) a21=0.5d0*chemof1(DKbs,S1(i,j))*ub1(i,j) a21=a21+0.5d0*chemof1(DKbs,S1(1—1,j))*ub1(i-1,j) b12=0.5d0*chemof1(DKbs,S1(1,j+1))*ub1(1,j+1) b12=b12+0.5d0*chemof1(DKbs,S1(i,j))*ub1(i,j) b21=0.5d0*chemof 1(DKbs,S 1(i,j))*ub1(i,j) b21=b21+0.5d0*chemof1(DKbs,S 1(1,j- 1))*ub1(i,j— 1) c12=0.5d0*chemof1(DKbq,ql(1+1,j))*ub1(i+1,j) ’ c12=c12+0.5d0*chemof1(DKbq,q1(i,j))*ub1(i,j) 021:0.5d0*chemof1(DKbq,q1(i,j))*ub1(i,j) c21:021+0.5d0*chemof1(DKbq,q1(1~1,j))*ub1(i-1,j) d12=0.5d0*chemof1(DKbq,q1(i,j+1))*ub1(i,j+1) d12=d12+0.5d0*chemof1(DKbq,q1(i,j))*ub1(i,j) d21=0.5d0*chemof1(DKbq,q1(i,j))*ub1(i,j) d21=d2 1+0.5d0*chemof1(DKbq,q1(i,j-1))*ub1(i,j- 1) rhsub(i)=ub1(i,j)+D6*(ub1(i+1,j)-ub1(i,j)) * +D6*(ub1(i-1,j)-ub1(i,j)) rhsub(i)=rhsub(j)+(ht/2.0d0)*f(cb,vb,h1(i,j))*ub1(i,j) rhsub(i)=rhsub(i)-D7*a12*(s1(1+1,j)-sl(1,j)) rhsub(j)=rhsub(j)-D7*a21*(sl(i—1,j)-s1(i,j)) rhsub(i)=rhsub(i)-D7*b12*(s1(i,j+1)-sl(i,j)) rhsub(i)=rhsub(j)-D7 *b2 1 *(s 1(1,j- 1)—s1(i,j)) rhsub(i)=rhsub(j)-D8*c12*(q1(i+1,j)-q1(i,j)) rhsub(i)=rhsub(i)-D8*021*(ql(1-1,j)-q1(1,j)) rhsub(i)=rhsub(j)-D8*dl 2*(q1(i,j+1)-q1(i,j)) rhsub(i)=rhsub(j)~D8*d21*(q1(i,j-1)-q1(i,j)) funub(j)=1.0d0 1080 continue C call PDMTRIX(m,D 1,funs,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhss,us) call PDMTRIX(m,D1,funq,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhsq,uq) call PDMTRIX(m,D0,funh,aa,bb,cc) 178 call TRIS OLV(m,aa,bb,cc,rhsh,uh) call PDMTRD{(m,D2,funua,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhsua,uua) call PDMTRIX(m,D6,funub,aa,bb,cc) call TRISOLV(m,aa,bb,cc,rhsub,uub) do 1085 j=0,m 32(i,j)=usti) q2(i.j)=uqti) h2(i,j)=uh(i) ua2(i,j)=uua(i) ub2(i,j)=uub(i) 1085 continue 1090 continue C do 1095 j=0,m do 1096 i=0,m 81(i.j)=sz q1(i,j)=q2(i,j) h1(i,j)=h2(i,j) ua1(i,j)=ua2(i,j) ub1(1,j)=ub2(i,j) 1096 continue 1095 confinue C 2000 continue C call veloc(mvel,lvel,distm,t1mm,distl,t1ml,velm,ve11) write(30,*)'timemax,distmax,velm' write(3 1,*)'timelim,distlim,vell' kvel2=mcount2lnumvel do 2060 i=1,mvel—1 write(30,*)timm(i),',',distm(i) ,',',ve1m(1) 2060 continue do 2062 1: 1 ,lvel-l write(31,*)timl(i),',',distl(i),',',vell(1) 2062 continue if(movc.eq.2)then write(16,*)'];' write(l7,*)'];’ write(23,*)'] ;’ write(18,*)'];' endif stop end 179 C C C The following four functions give the initial conditions C C function ua0(x,y,ua1ni,yoff,xoff,aw1d) implicit double precision (a-h,o-z) 301 format(6x,d10.9) r=sqrt((x-xoff)**2+(y—yoff)**2) if(r.1e.0.3d0) then ua0=uaini/dexp(awid*r) else ua0=0.0d0 endif return end function ub0(x,y,ubini,yoff,xoff,bwid) implicit double precision (a—h,o-z) 301 format(6x,d10.9) r=sqrt((x-xoff)**2+(y-yoff)**2) if(r.le.0.3d0) then ubO=ub1ni/dexp(bwid*r) else ub0=0.0d0 endif return end C C C The function f(x) is the uptake of nutrient C function f(ca,v,x) implicit double precision (a-h,o-z) if(x.ge.0.0d0) then f=v*x/(ca+x) else f=0.0d0 endif return end C C fs=f(x)/x function fs(c,v,x) 180 implicit double precision (a-h,o-z) if(x.ge.0.0d0) then fs=v/(c+x) else fs=0.0d0 endif return end C C chemof1*chemof2: non-linear chemotaxis coefficient function C chemofl is defined by Keller-Segel model C function chemof1(DK,x) implicit double precision (a-h,o-z) if(x.ge.0.0d0) then chemof1=DK/(DK+x)**2 else .chemof1=0.0d0 endif return end C C function diffu is defined by Monod's law C function diffu(v2,c2,eps,x) implicit double precision (a—h,o-z) if(x.ge.0.0d0) then diffu=v2*x/(02+x)+eps else diffu=eps endif return end C tr(i): sets tr to 0 if outside window, or to desired value inside window function tr(i,h,tk,dxy) implicit double precision (a—h,o-z) x=-2.5d0+1*h y=abs(x) if(y.le.2.0d0) then tr=tk/dxy else tr=0.0d0 endif return 181 end C C rhbound: subroutine to calculate boundary conditions at C sink or reservoir. C subroutine rhbound(n,a,b,c,d,hr,t,rh) implicit double precision (a-h,o—z) C dimension a(0:n),b(0:n),rh(0:n) rh=a+d*(b-a-hr*t*(a-c)) return end C rh2: subroutine to calculate boundary conditions at C sink or reservoir. C subroutine rh2(n,a,b,e,c,d,hr,t,rh) implicit double precision (a-h,o-z) C dimension a(0:n),b(0:n),rh(0:n) rh=a+d*(b-2*a+e+hr*t*c) return end C C maxcalc: finds maximum density on wave C subroutine maxcalc(ua1,k,m,ht,hr,umax,pmax,t) implicit double precision (a—h,o-z) dimension ua1(0: 100,0: 100) mcent=m/2 1max=0 5 if(ua1(mcent,imax+1).ge.ua1(mcent,imax))then imax=imax+1 goto 5 endif umax=ua1 (mcent,imax) pmax=2.5d0'-float(imax)*hr C write(*,*)'pmax=',pmax =k*ht C write(*,*)'t=',t return end C C limcalc: finds first point where minimum cell density occurs C subroutine limcalc(ua1,k,m,ht,hr,dmin,ulim,plim,t) implicit double precision (a-h,o-z) ' 182 dimension ua1(0: 100,0: 100) mcent=m/2 imax=0 5 if(ua1(mcent,imax+1).1t.dmin)then imax=imax+1 goto 5 endif ulim=ua1(mcent,imax) plim=2.5d0-float(imax)*hr C write(*,*)'plim=',plim t=k*ht C write(*,*)'t=',t return end subroutine veloc(mvel,lve1,distm,timm,distl,timl,velm,vell) implicit double precision (a—h,o-z) dimension distm(0: 100),distl(0: 100) dimension t1mm(0: 100),timl(0: 100) dimension velm(0: 100),ve11(0: 100) do 2089 i=2,mvel-1 velm(i—1)=(distm(i)-distm(i-1))/(timm(i)-timm(i- 1)) 2089 continue do 2099 i=2,lvel-1 vell(i-1)=(distl(i)-distl(i-1))/(timl(i)-timl(i- 1)) 2099 continue return end C C SOLVES THE TRIDIAGONAL SYSTEM OF LINEAR EQUATIONS C SUBROUTINE TRIS OLV(N,A,B,C,FUN,X) implicit double precision (a—h,o-z) DIMENSION A(O:N),B(0:N),C(O:N),FUN(0:N),X(0:N) DIMENSION ARFA(0:500),BATA(0:500),GAMA(O:500),GUN(0:500) B(O)=0.0d0 C(N)=0.0d0 ARFA(O)=A(0) GAMA(O)=C(0) DO 4000 I=1,N GAMA(I)=C(I) BATA(D=B (D/ARFA(I- 1) 183 ARFA(I)=A(I)-BATA(I) *C(I— 1) 4000 CONTINUE GUN(O)=FUN(0) DO 4999 I=1,N GUN(I)=FUN(I)-BATA(I)*GUN(I- 1) 4999 CONTINUE X(N)=GUN(N)/ARFA(N) DO 4888 I=1,N X(N-I)=(GUN(N-I)—GAMA(N-I)*X(N-I+1))/ARFA(N-I) 4888 CONTINUE C RETURN END C C SET UP THE TRIDIAGONAL MATRD( for constant coefficient terms C SUBROUTINE PDMTRIX(n,dxy,fu,aa,bb,cc) implicit double precision (a-h,o-z) DINIENSION fu(0:n),aa(0:n),bb(0:n),cc(0:n) bb(0)=0.0 aa(0)=fu(0)+1 .0d0*dxy cc(0)=-1 .0d0*dxy bb(n)=-1 .0d0*dxy aa(n)=fu(n)+1.0d0*dxy cc(n)=0.0 do 5000 i=1,n-1 bb(i)=-dxy cc(i)=-dxy aa(i)=fu(i)+2.0d0*dxy 5000 continue return end APPENDIX E APPENDIX E Input files for model the .. file 5615 I! am . 1 E. The Par ' Sample in put files for the program given 1n Appendlx f1 5615 the initial . Il'nit" ] 6 values of all parameters and the time for solution. The 1 cifiCS which , ll ' 6 conditions for each population and chemical component The tfine—points should be printed to the output files. Sample param file: In: 41 t1minit= 06.0d0 timin0c= 30.0d0 numcel: 2 numvel: 100 kstep: 20 dmin = 0.000001d0 . move: 1 1=individual matrices, 2=m0V15 R =. 2.5d0 ayoff: 0.0d0 axoff: 0.0d0 byoff= 0.0d0 bxoff: 0.0d0 awid = 1.0d0 bwid = 1.0d0 Dh — 0.01d0 thn 0.155d0 ths 0.155d0 the 0.000d0 thw 0.000d0 Ds 0.033d0 tsn 0.1071d0 tss 0.1071d0 tse 0.000d0 tsw 0.000d0 Dacs = 0.03d0 Dbcs = 0.00d0 DKas = 0.00000200d0 184 185 DKbs = 0.00000200d0 cas = 0.0000550d0 vas = 0.600d0 Yas = 1.0000d0 cbs = 0.0000550d0 vbs = 0.000d0 Ybs = 1.0000d0 Dq = 0.033d0 tqn = 0.891d0 tqs = 0.891d0 tqe = 0.000d0 tqw = 0.000d0 Dacq= 0.080d0 Dbcq = 0.00d0 DKaq = 0.0000330d0 DKbq = 0.0000330d0 caq = 0.000000067d0 vaq = 0.02d0 Yaq = 1.0000d0 cbq = 0.000000067d0 qu = 0.0d0 qu = 1.0000d0 Dua = 0.00100d0 va = 0.35d0 ca = 0.00000408d0 Ya = 0.350d0 Dub = 0.0010d0 vb = 0.50d0 cb = 0.00000408d0 Yb = 0.50d0 Sample init file: uaini: 0.000005d0 ubini= 0.000005 d0 rsn = 0.000132d0 rss = 0.00000d0 rse = 0.0000d0 rsw = 0.0000d0 rqn = 0.0000132d0 rqs = 0.0000132d0 rqe = 0.0000132d0 rqw = 0.0000132d0 186 rhn = 0.00046d0 rhs = 0.00046d0 the = 0.00046d0 rhw = 0.00046d0 stini= 0.00000d0 qtini: 0.0000132d0 htini= 0.00046d0 Sample time file: numti=04 time1= 10.0d0 time2= 15.0d0 time3= 20.0d0 time4= 30.0d0 APPENDIX F APPENDIX F Matlab program files The following Matlab program files were created for use with Matlab 4.2c for Windows. No guarantees are made that they will work with more recent versions. Commands enclosed in quotes should be entered into Matlab exactly as typed. The file appmultm plots the cell matrix as a density plot. The simulation plots in Figure 13 were created using appmultm, for example. The syntax is "appmult(u,#,'pcolor')" where u is the name of the matrix to be plotted, and # is the amount of interpolation; typically 2 or 3 works well. (Note: I did not write this file. It is modified from the version found at the Mathworks site, www.mathworks.com.) function hfi = app_int(x,y,z,s,dc) % This function approximates interpolated shading by interpolating % data and using flat shading. The inputs are the data that was used % to generate the original object, a scale factor, and a string % that contains the drawing command. The % function returns a handle to the new % object. % % Syntax 1:For just Z—Data (e.g., surf(z) => app_int(z,s,dc)) % z = zdata, s = scaling factor and dc = drawing command used to % create the plot % % Syntax 2: For x,y,z Data ( e.g., surf(x,y,z) ) => app_int % (x,y,z,s,dc) % % Example 1: [x,y,z] = peaks;surf(z);shading interp % app_int(z,3,'surf) interpolate by a factor of 3 % % Example 2: [x,y,z ] = peaks; surf(x,y,z); shading interp % app_int(x,y,z,3,‘surf')interpolate by a factor of 3 if ( nargin == 3 ) do = z; 187 188 z = x; s = y; [m n] = size(z); x = 1:n; y = (1:m)'; mscal = m*s; nscal = n*s; xi = linspace(1,n,nscal); yi = linspace(1,m,mscal)'; 21 = interp2(x,y,z,xi,yi); % figure; cmd = ['hfi=' dc '(xi,yi,zi);']; i eval(cmd); shading flat; elseif ( nargin == 5 ) [m n] = size(z); mscal = m*s; nscal = n*s; mint = min(min(x)); maxt = max(max(x)); xi = linspace(mint,maxt,mscal); mint = min(min(y)); maxt = max(max(y)); yi = linspace(mint,maxt,nscal)'; zi = interp2(x,y,z,xi,yi); figure; cmd = ['hfi=' dc '(xi,yi,zi);']; eval(cmd); shading flat; end colormap(gray) axis('square') The program surftwom allows two cell matrices to be plotted on the same figure, with different colormaps. Examples of this program are Figure 33C and Figure 33D. The syntax is "surftwo(u1,u2)" where 111 and u2 are the two cell matrices. %Program to plot two cell populations on one figure %with different colormaps. function[h]=surftwo(ua,ub) h(1)=surf(ua),axis equal 189 hold on h(2)=surf(ub); hold off m=64; cmin=min(ua(: )); cmax=max(ua(: )); cl=min(m, round((m-1)*(ua- -cmin)/(cmax- -cmin))+1); cmin=min(ub(: )); cmax=max(ub(:)); c2=min(m,round((m-1)*(ub-cmin)/(cmax-cmin))+1)+64; set(h(1),'cdata',cl); set(h(2),'cdata',c2); caxis([min(cl(:)) maX(c2(:))]); view(2); shading interp; colormap([bone(64);copper(64)]); The file centmassm calculates the average mass for two cell matrices. An example is Figure 35F. The syntax is "[cma,cmb]=centmass(u1,u2)" where 111 and u2 are the cell matrices. The returned variables cma and cmb may be plotted with "plot(cma), hold, plot(cmb)". function[cma, cmb]=centmass(ua, ub) %calculate the total mass along line perpendicular to gradient, %by summing across each line parallel to gradient. [r,c]=size(ua); centline=round(r/2); for 1: 1 :r; cma(i)=sum(ua(i,:))/(r- 1); cmb(i)=sum(ub(i,:))/(r—1); end ' The file competit.m computes the dynamic competition factor for two cell populations. Examples are Figure 35A and Figure 35B. The syntax is "competit", where 190 the value for num (the number of timepoints) and the cell matrices, have already been entered into the Matlab variable space. The returned variable rab is the dynamic competition factor, and is the ratio of the total mass of Population A to the total mass of Population B. A plot can be created using "plot(rab)". %Computes total mass of a and b and ratio a/b for uai, ubi where i=1:num %To use, first specify value for num, then type competit clear ta,clear tb,clear rab for i=1:num; [ta(i)]=totmass(eva1( ['ua',num28tr(i)]), 1 .6); [tb(i)]=totmass(eval(['ub’,num28tr(i)]),1.6); rab(i)=ta(i)/tb(i); end The file totmass.m is a subroutine called by the program competit.m. function[tmass]=totmass(u,ht) % Calculates total mass under cell profile % Assumes length and width of DGC = 5 cm 4 % ht is the height of the gel, usually 1.6 cm % Find number of increments [r,c]=size(u); % Compute area under 2-D plot, multiply by l*w*h for mass tmass=trapz(trapz(u))*5/(r-1)*5/(c-1)*ht; The program totvelm calculates the mass fluxes due to chemotaxis to S, chemotaxis to Q, random motility, and the total mass flux. Examples are Figure 21-Figure 24. The syntax is "[uxs,uys,uxq,uyq,uxu,uyu,uxt,uyt,mag,dir,dxc,dyc,aoverh]=totvel(kds,xos,kdq,xoq,du,s, q,u,minv)" where kds and kdq are the dissociation constants for the receptor—S and 191 receptor-Q complexes, respectively; xos and xoq are the chemotactic sensitivity coefficients for S and Q; du is the random motility coefficient; 5, q, and u are the S- chemoattractant, Q-chemoattractant, and cell matrices, respectively; and minv is a minimum flux below which no values are reported. The returned variables uxs and uys are the x and y components of the flux due to chemotaxis to S; uxq and uyq are the x and y components of the flux due to chemotaxis to Q; uxu and uyu are the x and y components of the flux due to random motility; and uxt and uyt are the x and y components of the total mass flux. The fluxes can be plotted with the quiver command, as illustrated near the end of the program. function[uxs,uys,uxq,uyq,uxu,uyu,uxt,uyt,mag,dir,dxc,dyc,aoverh]=totvel(kds,xos ,kdq,xoq,du,s,q,u,minv) %This program calculates the flux matrices for chemotaxis to S and Q, %(uxs,uys,uxq,uyq), random motility (uxu,uyu), and the total flux (uxt,uyt). [m,n]=size(u); dir=zeros(m,n); %Calls to the subroutine chemo, which is where intermediate fluxes actually calculated. [uxs,uys]=chemo(kds,xos,s,u,minv); [uxq,uyq]=chemo(kdq,xoq,q,u,minv); [uxu,uyu]=randmot(du,u,minv); %Total flux found by summing components of flux. uxt=uxs+uxq+uxu; uyt=uys+uyq+uyu; %Magnitude of maximum flux calculated. mag=sqrt(uxt.*uxt+uyt.*uyt); %Maximum flux magnitude found. maxvel=max(max(abs(mag))); %A minimum flux value to be graphed is calculated. Minv=5 works well. (80% of fluxes graphed). minvel=maxvel/minv; %Loop to remove fluxes whose magnitudes are lower than minvel. fori = lzm, for j = 1:n; if abs(mag(i,j)) < minvel uxt(i,j)=0; 191 receptor-Q complexes, respectively; xos and xoq are the chemotactic sensitivity coefficients for S and Q; du is the random motility coefficient; 3, q, and u are the S- chemoattractant, Q-chemoattractant, and cell matrices, respectively; and minv is a minimum flux below which no values are reported. The returned variables uxs and uys are the x and y components of the flux due to chemotaxis to S; uxq and uyq are the x and y components of the flux due to chemotaxis to Q; uxu and uyu are the x and y components of the flux due to random motility; and uxt and uyt are the x and y components of the total mass flux. The fluxes can be plotted with the quiver command, as illustrated near the end of the program. function[uxs,uys,uxq,uyq,uxu,uyu,uxt,uyt,mag,dir,dxc,dyc,aoverh]=totvel(kds,xos ,kdq,xoq,du,s,q,u,minv) %This program calculates the flux matrices for chemotaxis to S and Q, %(uxs,uys,uxq,uyq), random motility (uxu,uyu), and the total flux (uxt,uyt). [m,n]=size(u); dir=zeros(m,n); %Calls to the subroutine chemo, which is where intermediate fluxes actually calculated. [uxs,uys]=chemo(kds,xos,s,u,minv); [uxq,uyq]=chemo(kdq,xoq,q,u,minv); [uxu,uyu]=randmot(du,u,minv); %Total flux found by summing components of flux. uxt=uxs+uxq+uxu; uyt=uys+uyq+uyu; %Magnitude of maximum flux calculated. mag=sqrt(uxt.*uxt+uyt.*uyt); %Maximum flux magnitude found. maxvel=max(max(abs(mag))); %A minimum flux value to be graphed is calculated. Minv=5 works well. (80% of fluxes graphed). minvel=maxvel/minv; %Loop to remove fluxes whose magnitudes are lower than minvel. for 1 = l:m, for j = 1:n; if abs(mag(i,j)) < minvel uxt(i,j)=0; 192 uyt(i,j)=0; end end end %subplot(2,2,1),quiver(uxs,uys,3,'w'),title('chemo s'),axis square %subplot(2,2,2),quiver(uxq,uyq,3,'w'),title('chemo q'),axis square %subplot(2,2,3),quiver(uxu,uyu,3,'w'),t1t1e('randmot'),axis square %subplot(2,2,4),quiver(uxt,uyt,3,‘w'),title('total flux'),axis square dxc=zeros(m,n); dyc=zeros(m,n); for i = l:m, for j = 1:n; if mag(i,j)~=0 dxc(i,j)=uxt(i,j)/mag(i,j); dycti,j)=uyt(i,j)/mag(i,j); end end end The file chemo.m is called as a subroutine 1n totvel.m. function[ux,uy,satx,saty]=chemo(kd,xo,c,u,minv) %calculate the partial derivatives of c in the x and y directions [m,n]=size(c); dx=5/m; dy=5/n; [px,py]=gradient(c,dx,dy); %calculate the matrix kd/(kd+c)"2 for i = l:m, for j = 1:n; sat(i,j)=kd/(kd+c(i,j))"2; end end %calculate the velocity coefficient matrix kd/(kd+c)"2 *px(or py) %the .* operator causes only the corresponding two matrix positions to be multiplied %rather than real matrix multiplication satx=xo*(sat.*px); saty=xo*(sat.*py); %clear all values of u below a minimum value maxu=max(max(u)); 193 minu=maxul 1000; for i = l:m, for j = 1:n; if u(i,j) < minu u(i,j)=0; end end end %multiplies the cell density times the velocity coefficient ux=u.*satx; uy=u.*saty; mag=sqrt(ux.*ux+uy.*uy); maxvel=max(max(abs(mag))); minvel=maxve1/minv; %maxvely=max(max(abs(uy))); %minvely=maxvely/minv; fori = l:m, for j = 1:n; if abs(mag(i,j)) < minvel ux(i,j)=0; uy(i,j)=0; end end end The file randmotm is called as a subroutine in the program totve1.m function[ux,uy]=randmot(du,u,minv) %calculate the partial derivatives of u in the x and y directions [px,py]=gradient(u); %calculate the velocity coefficient matrix du*px(or py) %the .* operator causes only the corresponding two matrix positions to be multiplied %rather than real matrix multiplication %clear all values of 11 below a minimum value [m,n]=size(u); maxu=max(max(u)); minu=maxul 1000; for i = l:m, for j = 1:n; if u(i,j) < minu 194 u(i,j)=0; end end end ux=du*u.*px; uy=du*u.*py; mag=sqrt(ux.*ux+uy.*uy); maxvel=max(max(abs(mag))); minvel=maxvel/minv; fori = l:m, for j = 1:n; if abs(mag(i,j)) < minvel ux(i,j)=0; UY(i,j)=O; end end end The program "retard4.m" calculates the global chemotactic response factor. An example is Figure 25. The syntax is "[phi,den,num,dent]=retard4(u,s,kd,perct)" where u is the cell matrix, s is the chemoattractant S matrix, kd is the dissociation constant for the receptor-S complex, and perct is the lowest percent value of the maximum cell concentration for which a value of the global chemotactic response factor will be calculated. The global chemotactic response factor is given by the returned variable phi. function[phi,den,num,dent]=retard4(u,s,kd,perct) % Calculation of retardation factor. % Retard4 calculates factor along y-direction for all x. % perct is the minimum percentage of the maximum cell concentration above % which the cell concentration will not be assumed to be 0. % Right now, only works when gradient in direction of matrix columns. % Calculation of s-gradient [m,n]=size(u); dy=5/m; maxui=max(max(u)); minui=perct*maxui; 195 for i=1:n, grads=abs(gradient(s(:,i)',dy)'); % Calculate (1+s/kd)) term. grp1=(1+s(:,i)/kd); % Calculate (1+s/kd)"(-2)*gradient(s) gsgrp1=grp1."(-2).*grads; % Calculate integrand of numerator of retardation factor num=gsgrp1.*u(:,i); for j: 1 :m, if u(i,i) < minui num(j)=0; end end % Calculate integrand of denominator of retardation factor den=grads.*u(:,i); % Evaluate integrals with trapezoid method. Calculate phi. % phi=zeros(size(pos)); %Make den a function of i, so average maximum flux can be reported. dent(i)=trapz(den); phi(i)=trapz(num)/trapz(den); end APPENDIX G APPENDIX G Dimensionless model Over the course of the competition modeling studies, it became apparent that there was an almost infinite number of parameter combinations to be tested. Many times, it is possible to reduce the number of parameters in a model by making it dimensionless, and then varying the dimensionless groups. The competition model was made dimensionless by employing the following scales (the ' indicates the dimensional quantity): t = vnHt’ (57) -' 58 ui = u, for 1=a,b ( ) uiO . °’ . (59) _]=—_" forj=H,S,Q .10 D (60) V = V’ —-H- VaH When these scales were substituted into the dimensional equations and the equations were rearranged, the new dimensionless equations were: 196 197 a (61) But“ :1, Vzua—X35$V.[[___G_IS_2]URVS —}ea83QV. _Q‘Q—z. 11an + H 11,, (Gas +8) (eaQ + Q) Wait '1' H Bub 0 9 H (62) =?\. Vzu —?»5 V- —LuVS —?\.5 V— ——bQ—uV + n" at b b b w [[(o,,+s)2] b b bQ (9,,Q+Q)2 b Q WbH+Hub as S S (63) —=o V28— at S Bus Was-1's ua_ 1tl-IBbS—I—Wbs"1'SUb 519_ 2 __9_ Q (64) at —GQV Q BaQ WaQ'l'Qua—TCHBbQ wa+Qub 8H H H (65) __ ____ V2 __ _ __— at H 1’ tllafl+Hua RH¢b tllbfl+Hub where the dimensionless variables are defined as: - . . . . . 66 l, = E‘— for i=a,b Measure of random motllity to nutr1entd1ffus1on ( ) H Xoij . . . . . . _ ( 67) 5,1. = —- for 1=a,b; j=S,Q Measure of chemotaXIS to random 1not1t111ty Hij .. ( 68) 0,]. - ,le for i=a,b Jo 198 \Ilij = "fl for i=a,b; j=S,Q Saturation constant over initial concentration 0 V . 71:]- = —b’ for j=S,Q,e,d Max. growth rate of b on j to max. growth rate of a on j at Vi'uio . . , . . 13,,- = J— for 1=a,b; j=S,Q,H,e Compares max uptake of1 on j to rate of 1 on H 1H 0 .l D. o. = D—’ for j=S,Q Ratio of diffusivity of j to diffusivity of H H (l)i = '0 for 1=a,b; St01chiometr1c concentration ratio HOYiI-l ( 69) ( 70) (71) ( 72) ( 73) By comparing the dimensionless equations (Equations ( 61)-( 65)) to the dimensional equations (Equations (38)-(40)), it is obvious that the number of parameters has been reduced by only one. The introduction of the dimensionless groups has not significantly simplified the model. 198 ij . . . . . . . “’0' = ,— for i=a,b; j=S,Q Saturation constant over initial concentration 0 v . it]. = —b’ for j=S,Q,e,d Max. growth rate of b on j to max. growth rate of a on j at i'uio . . [3,). = J for i=a,b; j=S,Q,H,e Compares max uptake of1 on j to rate of1 on H 1H0 l D. G. = D-J- for j=S,Q Ratio of diffusivity of j to diffusivity of H H d), = ‘0 for 1=a,b; Steichiometric concentration ratio HOYiH ( 69) ( 70) (71) ( 72) ( 73) By comparing the dimensionless equations (Equations ( 61)-( 65)) to the dimensional equations (Equations (38)—(40)), it is obvious that the number of parameters has been reduced by only one. The introduction of the dimensionless groups has not significantly simplified the model. LIST OF REFERENCES LIST OF REFERENCES Abu-Ashour, J., D. M. Joy, H. Lee, H. R. Whiteley, and S. Zelin. 1994. Transport of Microorganisms through Soil. Water, Air and Soil Pollution. 75: 141—158. Adler, J. 1966. Chemotaxis in Bacteria. Science. 153, 708—716. Adler, J. 1972. A Method for Measuring Chemotaxis and Use of the Method to Determine Optimum Conditions for Chemotaxis by Escherichia coli. Journal of General Microbiology. 74: 77-91. Agladze, K., L. Budriene, G. Ivanitsky, V. Krinsky, V. Shakhbazyan, and M. Tsyganov. Wave mechanisms of pattern formation in microbial populations. Proc. R. Soc. Lond. 253: 131—135. Ames, T. 1997. A Continuous Magnetofluidized Bed Bioreactor for Production of Plant— Cell Secondary Metabolites. Doctoral Dissertation. Michigan State University. Bailey, 1. E. and D. F. 0111s. 1986. Biochemical Engineering Fundamentals. 2'“l edition. McGraw—Hill Inc., New York. Barton, J. W. and R. M. Ford. 1997. Mathematical Model for Characterization of Bacterial Migration through Sand Cores. Biotechnology and Bioengineering. 53: 487 -496. Berg, H. C. 1988. A Physicist Looks at Bacterial Chemotaxis, p. 1-9. In Cold Spring Harbor Symposia on Quantitative Biology, Vol. III. Berg, H. C. and P. M. Tedesco. 1975. Transient response to chemotactic stimuli in Escherichia coli. Proc. Nat. Acad. Sci. USA. 72: 3235-3239. Berg, H.C. and D. A. Brown. 1972. Chemotaxis in Escherichia coli Analysed by Three- Dimensional Tracking. Nature. 239:500—504. Boon, 1., and B. Herpigny. 1986. Model for chemotactic bacteria bands. Bull. Math. Biol 48: 1—19. Bosma, T. N. P., J. L. Schnoor, G. Schraa, and A. J. B. Zehnder. 1988. Simulation Model for Biotransformation of Xenobiotics and Chemotaxis in Soil Columns. J. of Contaminant Hydrogeology. 2: 225-236. Budrene, E. O. and H. C. Berg. 1995. Dynamics of formation of symmetrical patterns by chemotactic bacteria. Nature. 376: 49—53. Caldwell, D. E. and P. Hirsch. 1972. Growth of microorganisms in two—dimensional steady—state diffusion gradients. Can. J. Microbiol. 19: 53—58. Carnahan, B., H. A. Luther, and I . O. Wilkes. Approximation of the Solution of Partial Differential Equations, p. 429-464. In Applied Numerical Methods. John Wiley and Sons, Inc. 1969. 199 200 Chapra, S. C. and R. P. Canale. Parabolic Equations in Two Spatial Dimensions, p. 745— 748. In Numerical Methods for Engineers. McGraw-Hill Book Company. 1988. Chet, I. and R. Mitchell. 1976. Ecological Aspects of Microbial Chemotactic Behavior. Ann. Rev. Microbiol. 302221-39. Criddle, C. S., J. T. DeWitt, D. Grbic-Galic, and P. L. McCarty. 1990. Transformation of Carbon Tetrachloride by Pseudomonas sp. Strain KC under Denitrification Conditions. Appl. Environ. Microbiol. 56:3240—3246. Cussler, E. L. 1984. Diffusion: Mass transfer in fluid systems. 1St edition. Cambridge University Press, Cambridge. Devare , M. and M. Alexander. 1995. Bacterial Transport and Phenanthrene Biodegradation in Soil and Aquifer Sand. Soil. Sci. Am. J. 59: 1316—1320. Duffy, K. J., P. T. Cummings, and R. M. Ford. 1995. Random Walk Calculations for Bacterial Migration in Porous Media. Biophysical Journal. 68:800—806. Dybas, M.J., G.M. Tatara, and CS. Criddle. 1995. Localization and Characterization of the Carbon Tetrachloride Transforming of Pseudomonas sp. strain KC. Appl. Environ. Microbiol. 61:758—762. Emerson, D., R. M. Worden, and J. A. Breznak. 1994. A Diffusion Gradient Chamber for Studying Microbial Behavior and Separating Microorganisms. Appl. Environ. Microbiol. 60,1269-1278. Emerson, D., S. F. Peteu, and R. M. Worden. 1996a. A Catalase Microbiosensor for Detecting Hydrogen Peroxide. Biotechnology Techniques. 10:673-678. Emerson, D., S. Peteu, M. Widman, and M. Worden. 1996b. A Microbiosensor for Sensing Glucose, Galactose, or Choline, and its Application for Measuring Glucose Gradients in Semi—solid Gels. The Sixth International Meeting on Chemical Sensors. NIST, Gaithersburg, MD. July 22-25, 1996. Ford, R. M. 1992. Mathematical Modeling and Quantitative Characterization of Bacterial Motility and Chemotaxis, p. 177-215. In C. J. Hurst (ed), Modeling the Metabolic and Physiologic Activities, of Microorganisms. John Wiley & Sons, Inc. New York. Ford, R. M. and D. A. Lauffenburger. 1990. Measurement of Bacterial Random Motility and Chemotaxis Coefficients: H. Application of Single—Cell—Based Mathematical Model. Biotech. Bioeng. 37: 661-672. Ford, R. M., B. R. Phillips, J. A. Quinn, and D. A. Lauffenburger. 1990. Measurement of Bacterial Random Motility and Chemotaxis Coefficients: l. Stopped-Flow Diffusion Chamber Assay. Biotech. Bioeng. 37: 647—660. Frymier, P. D., R. M. Ford and P. T. Cummings. 1992. Cellular Dynamics Simulations of Bacterial Chemotaxis. Chemical Engineering Science. 48:687-699. 201 Frymier, P. D., R. M. Ford and P. T. Cummings. 1994. Analysis of Bacterial Migration: 1. Numerical Solution of Balance Equation. AIChe Journal. 40: 704-715. Hansen, S. R. and S. P. Hubbell. 1980. Single-Nutrient Microbial Competition: Qualitative Agreement Between Experimental and Theoretically Forecast Outcomes. Science. 207: 149 1- 1493. Harwood, C. S., R. E. Porales, and M. Dispensa. 1990. Chemotaxis of Pseudomonas putida toward Chlorinated Benzoates. Appl. Environ. Microbiol. 56: 1501—1503. Holmes, E. E., M. A. Lewis, J. E. Banks, and R. R. Veit. Partial Differential Equations in Ecology: Spatial Interactions and Population Dynamics. Ecology. 75: 17-29. Kato, J., A. Ito, T. Nikata and H. Ohtake. 1992. Phosphate Taxis in Pseudomonas aeruginosa. J. Bacteriol. 174: 5149-5151. Keller, E. F. and L. A. Segel. 1971. Travelling Bands of Chemotactic Bacteria: a Theoretical Analysis. J. Theor. Biol. 30: 235—248. Kelly, F. X., K. J. Dapsis, and D. A. Lauffenburger. 1988. Effect of Bacterial Chemotaxis on Dynamics of Microbial Competition. Microb. Ecol. 16: 115-131. Kennedy, M. J. and J. G. Lawless. 1985. Role of Chemotaxis in the Ecology of Denitrifiers. Appl. Environ. Microbiol. 49: 109-1 14. Lauffenburger, D. A., R. Aris, and K. Keller. 1982. Effects of Cell Motility and Chemotaxis on Microbial Population Growth. Biophys. J. 40: 209-219. Lauffenburger, D., and B. Calcagno P. 1983. Competition Between Two Microbial Populations in a Nonmixed Environment: Effect of Cell Random Motility. Biotech. and Bioeng. 25: 2130—2125. Lauffenburger, D., R. Aris, and K. H. Keller. 1981. Effects of random motility on growth of bacterial populations. Microb. Ecol. 7:207 -227 . Liu, Z. and K. D. Papadopoulos. 1995. Unidirectional Motility of Escherichia coli in Restricitve Capillaries. Applied and Environmental Microbiology. 61: 3567-3572. Liu, Z., W. Chen, and K. D. Papadopoulos. 1996. Bacterial Motility, Collisions, and Aggregation in a 3-1im-Diameter Capillary. Biotechnology and Bioengineering. 53: 14. Macnab, R. M. 1987. Motility and Chemotaxis, p. 732—759. In F. C. Neidhardt, J. L. Ingraham, K. B. Low, B. Magasanik, M. Schaechter, and H. E. Umbarger (ed), Escherichia coli and Salmonella typhimurium: Cellular and Molecular Biology, vol. 1. American Society for Microbiology, Washington, D. C. Maniatis, T., E. F. Fritsch, and J. Sambrook. 1982. "Molecular Cloning: A Laboratory Manual." Cold Spring Harbor Laboratory, p. 68. 202 Mayotte, T.]., M.]. Dybas, and CS. Criddle. 1996. Bench-Scale Evaluation of Bioaugmentation to Remediate Carbon-Tetrachloride Contaminated Aquifer Materials. Ground Water, 34: 358-367. Mesibov, R., G. W. Ordal, and J. Adler. 1973. The Range of Attractant Concentrations for Bacterial Chemotaxis and the Threshold and Size of Response over This Range. J. Gen. Physiol. 62: 203-223. Mikola, Mark. 1996. Optimization of Biocatalytic 3-Dehydroshikimic Acid Production from D-Glucose in Escherichia coli. Master's Thesis. Michigan State University. Nikata, T., K. Sumida, J. Kato and H. Ohtake. 1992. Rapid Method for Analyzing Bacterial Behavioral Responses to Chemical Stimuli. Appl. Environ. Microbiol. 58: 2250-2254. Nossal, R. 1972. Boundary Movement of Chemotactic Bacterial Populations. Mathematical Biosciences. 13: 397-406. Perry, R. H. and D. Green (eds). 1984. Perry’s Chemical Engineers’ Handbook, 6th Edition. McGraw-I-Iill Book Company, pp. 3-258 to 3-287. Peteu, S. F., D. Emerson, and R. M. Worden. 1996. A Clark-type oxidase enzyme-based amperometric microbiosensor for sensing glucose, galactose, or choline. Biosensors and Bioelectronics. 11:1059-1071. Revsbech, N. P. 1989. An oxygen microsensor with a guard cathode. Limnol. Oceanogr. 34: 474—478. Rivero, M. A., R. T. Tranquillo, H. M. Buettner and D. A. Lauffenburger. 1989. Transport Models for Chemotactic Cell Populations Based on Individual Cell Behavior. Chemical Engineering Science. 44: 2881-2897. Rivero-Hudec, M. and D. A. Lauffenburger. 1986. Quantification of Bacterial Chemotaxis by Measurement of Model Parameters Using the Capillary Assay. Biotechnology and Bioengineering. 28: 1178-1190. Sarkar, A. K., G. Georgiou, and M. M. Sharma. 1994a. Transport of Bacteria in Porous Media: 1. An Experimental Investigation. Biotech. and Bioeng. 44: 489-497. Sarkar, A. K., G. Georgiou, and M. M. Sharma. 1994b. Transport of Bacteria in Porous Media: II. A Model for Convective Transport and Growth. Biotech. and Bioeng. 44: 499-508. Schmidt, S., M. T. Widman, and R. M. Worden. 1997. A Laser-Diffraction Capillary Assay for Measuring Microbial Random Motility. Biotech. Techniques. 11: 423- 426. Segel, L. A. 1977. A Theoretical Study of Receptor Mechanisms in Bacterial Chemotaxis. J. Appl. Math. 32: 653665. Sitaraman, R. S., S. H. Ibrahim, and N. R. Kuloor. 1963. J. Chem. Eng. Data. 8:198. 203 Soby, S. and K. Bergman. 1983. Motility and Chemotaxis of Rhizobium meliloti in Soil. Applied and Environmental Microbiology. 46: 995-998. Staffeld, P. 0., D. A. Lauffenburger, and J. A. Quinn. 1987. Mathematical Analysis of Cell Transport Phenomena: Bacterial Chemotaxis in the Capillary Assay. Chem. Eng. Comm. 58: 339-351. Strauss, I, P. D. Frymier, C. M. Hahn, and R. M. Ford. 1995. Analysis of Bacterial Migration: II. Studies with Multiple Attractant Gradients. AIChE Journal. 41: 402-414. Tan, Y., J. T. Gannon, P. Baveye, and M. Alexander. 1994. Transport of bacteria in an aquifer sand: Experiments and model simulations. Water Resources Research. 30: 3243-3252. Tatara, G. M., M. J. Dybas, and C. S. Criddle. 1993. Effects of Medium and Trace Metals on Kinetics of Carbon Tetrachloride Transformation by Pseudomonas sp. Strain KC. Appl. Environ. Microbiol. 49: 109-114. Tilman, D. 1994. Competition and Biodiversity in Spatially Structured Habitats. Ecology. 75: 2-16. Widman, M. T., D. Emerson, C. C. Chiu, and R. M. Worden. 1996. Modeling Microbial Chemotaxis in a Diffusion Gradient Chamber. Biotech. and Bioeng. 55: 191—205. Wilke, C. R. and P. Chang. 1955. Correlation of Diffusion Coefficients in Dilute Solutions. AIChE J. 1: 264-270. Witt, ME. 1994. Development of a Laboratory-Scale Model Aquifer System to Monitor a Carbon-Tetrachloride-Transforming Zone by Pseudomonas sp. strain KC, MS Thesis, Michigan State University, East Lansing, MI. Wolfe, A. J., and H. C. Berg. 1989. Migration of Bacteria in Semisolid Agar. Proc. Natl. Acad. Sci. USA 86: 6973—6977. Woodward, D. E., R. Tyson, M. R. Myerscough, J. D. Murray, E. O. Budrene, and H. C. Berg. Spatio-Temporal Patterns Generated by Salmonella typhimurium. Biophysical Journal. 68: 2181-2189. Yamamoto, K. and Y. Imae. 1993. Cloning and Characterization of the Salmonella typhimurium—specific chemoreceptor Tcp for taxis to citrate from phenol. Proc. Natl. Acad. Sci. USA. 90: 217-221. Zhulin, I. B., E. H. Roswell, M. S. Johnson, and B. L. Taylor. 1997. Glycerol Elicits Energy Taxis of Escherichia coli and Salmonella Iyphimurium. J. Bacteriology. 179: 3196—3201 . 71 111111111111 312930168623 111111111111 1