ABSTRACT Title of dissertation: MATHEMATICAL MODELS AND SIMULATIONS OF PHOTOTAXIS AND CANCER IMMUNE INTERACTIONS Amanda Rae Galante Doctor of Philosophy, 2012 Dissertation directed by: Doron Levy Department of Mathematics and Center for Scienti c Computation and Mathematical Modeling The macroscopic dynamics created by complex systems of microscopic cells can be better understood using mathematical modeling. In this dissertation, we study the dynamics of two di erent biological systems: phototaxis exhibited by cyanobacterium Synechocystis sp. and cancer immune interactions as a ected by the protein B7-H1. Synechocystis sp., a common unicellular freshwater cyanobacterium, has been used as a model organism to study phototaxis, i.e. motion in the direction of a light source. This microorganism displays a number of additional characteristics such as delayed motion, surface dependence, and a quasi-random motion, where cells move in a seemingly disordered fashion instead of in the direction of the light source. These unexplained motions are thought to be modulated by local interactions between cells. In this work, we formulate a model of local interactions between phototactic cells in order to study the structure of their quasi-random motion. We present a stochastic dynamic particle system modeling interacting phototactic cells. We ex- tend our model of local interactions to include global forcing due to light. We also add an activation process of cells as they become a ected by the presence of light. We study the parameter space of our model by deriving a system of ordinary dif- ferential equations that describe the dynamics of the system in one dimension. The simulations of our model are consistent with experimentally observed phototactic motion. The second part of this dissertation focuses on the surface protein B7-H1. This protein, also called PD-L1 and CD274, is found on carcinomas of the lung, ovary, colon and melanomas but not on most normal tissues. B7-H1 has been experimentally determined to be an anti-apoptotic receptor on cancer cells, where B7-H1-positive cancer cells have been shown to be immune resistant, and in vitro experiments and mouse models have shown that B7-H1-negative tumor cells are signi cantly more susceptible to being repressed by the immune system. We derive a mathematical model for studying the interaction between cytotoxic T cells and tumor cells as a ected by B7-H1. By integrating experimental data into the model, we isolate the parameters that control the dynamics and obtain insights on the mechanisms that control apoptosis. MATHEMATICAL MODELS AND SIMULATIONS OF PHOTOTAXIS AND CANCER IMMUNE INTERACTIONS by Amanda Rae Galante Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial ful llment of the requirements for the degree of Doctor of Philosophy 2012 Advisory Committee: Professor Doron Levy, Chair Professor Radu Balan Professor Maria Cameron Professor Pierre-Emmanuel Jabin Professor Garegin Papoian, Dean?s Representative c Copyright by Amanda Rae Galante 2012 Acknowledgments During my time as a graduate student at the University of Maryland, I have been very fortunate to have been thoroughly supported academically, emotionally and nancially by many people. First, I must thank my academic advisor, Doron Levy. From the countless hours of critiquing my writing and providing priceless networking and travel oppor- tunities, to daily life and research advice, I thank you. I am very lucky to have had such a supportive, encouraging, and entertaining advisor. Many thanks to my dissertation committee, Radu Balan, Maria Cameron, Pierre-Emmanuel Jabin, and Garegin Papoian, for their careful reading and cri- tiquing of this work. Their insights and suggestions were most helpful. Other professors have provided immense guidance and support during my time in graduate school, particularly Eitan Tadmor and Konstantina Trivisa. Our collab- orators Devaki Bhaya, Susanne Wisen, and Koji Tamada provided many invaluable research discussions. I owe much gratitude to the many professors who instructed my graduate and undergraduate classes at the University of Maryland and Rose- Hulman Institute of Technology also helped to shape my academic career. I was nancially supported by the two-year Block grant fellowship, part-time teaching assistantships, and Doron?s research grants through NSF and NCI. I am also thankful for receiving the Ruth Davis fellowship and numerous travel grants that enabled my attendance at conferences, workshops, and summer schools. The friendship, advice, and numerous research-related and o -topic discus- ii sions with my colleagues, particularly, my academic twin Shelby Wilson, provided amusement, encouragement and academic support. I am grateful to have had awe- some, intelligent, friendly o cemates and academic associates, namely, Courtney Davis, Cristian Tomasetti, Prashant Athavale, and Anne Jorstad. Many thanks to James Greene and Daniel Weinberg for helping to carefully proofread this work. Many sta members of the Department of Mathematics and the Center for Scienti c Computation and Mathematical Modeling provided signi cant support. I truly appreciate the hard work of the CSCAMM sta , Agi Alipio, Je Henrikson, and Valerie Lum, who tolerated my countless o ce-related questions and requests. The support, friendship, and love of my family and friends made this expe- rience much more enjoyable. My friendships at the University of Maryland Space Systems Lab and in the Department of Mathematics helped me to adjust to grad- uate student life and East-coast living. I also enjoyed the refreshing company of a few Midwestern-transplant friends. Many hours of phone calls and emails with Katy Caruso helped to keep me sane. My parents, siblings, extended family and in-laws have provided unwavering emotional support. In particular, my parents have believed in me from the very beginning; the fruits of their labor have not gone unnoticed. My trouble-making pets, Zelda, Viktor, and Jekhar, provided daily dis- tractions and comfort. Finally, I thank my husband Joseph whose unconditional love, support, and toleration of many years of my craziness made this possible. My gratitude also goes to many unnamed others who helped me over the years. iii Table of Contents List of Tables vi List of Figures vii List of Abbreviations xvi 1 Introduction 1 1.1 Phototaxis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Biological background . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Mathematical models relevant to phototaxis . . . . . . . . . . 4 1.2 B7-H1 and tumor-immunology . . . . . . . . . . . . . . . . . . . . . . 5 1.2.1 Biological background . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Models relevant to B7-H1 and cancer immunology . . . . . . . 7 1.3 Making the connection between phototaxis and cancer . . . . . . . . 8 1.4 Overview of dissertation . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Modeling local interactions during motion of cyanobacteria 11 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Biological background . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 A mathematical model for local interactions of Synechocystis sp. . . . 20 2.3.1 Model assumptions . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.2 Model A: A local-interactions model . . . . . . . . . . . . . . 23 2.3.3 Model B: Random neighbor-dependent movement model . . . 25 2.4 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3 Phototaxis: incorporating global forcing due to the light source 40 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 Global forcing model for phototaxis . . . . . . . . . . . . . . . . . . . 42 3.2.1 Model assumptions . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2.2 Model H: Homogeneous response to directional light force . . . 45 3.2.3 Model I: Inhomogeneous response to directional light force . . 46 3.3 Simulation results of the models of Section 3.2 . . . . . . . . . . . . . 47 3.4 An activation model for phototaxis . . . . . . . . . . . . . . . . . . . 55 3.5 Simulation results of activation model . . . . . . . . . . . . . . . . . . 59 3.6 Conclusions and future directions . . . . . . . . . . . . . . . . . . . . 63 4 One-dimensional local interaction models 68 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2 Local interaction model for two populations in one dimension . . . . . 69 4.2.1 Derivation of the Reaction-Di usion Master Equation (RDME) 72 4.2.2 Deriving a system of ODEs . . . . . . . . . . . . . . . . . . . 79 4.3 An extended one-dimensional model of local interactions . . . . . . . 81 iv 4.3.1 Derivation of the Reaction-Di usion Master Equation . . . . . 82 4.3.2 Deriving a system of ODEs . . . . . . . . . . . . . . . . . . . 87 4.4 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4.1 Comparing the ODEs to the stochastic particle model . . . . . 91 4.4.2 Uniform initial conditions . . . . . . . . . . . . . . . . . . . . 97 4.4.3 Parameter analysis by simulation . . . . . . . . . . . . . . . . 101 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5 B7-H1 and Cancer Immune Interactions 118 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.2 A Mathematical model for CTL and tumor interaction . . . . . . . . 119 5.2.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.2.2 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.3 Expression for percent lysis . . . . . . . . . . . . . . . . . . . . . . . 125 5.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 6 General conclusions 134 Bibliography 137 v List of Tables 4.1 De nitions of particle-transposition operators acting on f~r;~l; ~rs; ~lsg. The rst subscript indicates the bin on which the operator is acting. The second subscript indicates the population type. The superscript indicates the direction in which that particle is moving, right (+), left (-), or (s) if the particle is becoming stationary. . . . . . . . . . . . . 85 5.1 Initial conditions for model . . . . . . . . . . . . . . . . . . . . . . . . 123 5.2 Parameter values for model simulation . . . . . . . . . . . . . . . . . 125 vi List of Figures 2.1 Experimental results illustrating (a) beginnings of nger-like pro- jection formation, (b) surface dependence, and (c/d) aggregation and nger formation. The bacteria within a matrix of extracellu- lar polysaccharides appear as white dots or circles. The bacteria are approximately 1 m in diameter. Each image contains wildtype cells exposed to a light source radiating from the upper left corner of each image. Inset (a) shows the front of a spot of cells 48 hours after plating. Inset (c) shows the exact same spot frame as (a), 24 hours later. Images (b) and (d) depict di erent sections of the same spot at di erent times. In (b), the density of cells is very low and surface de- pendence is observable|the cells are surrounded by polysaccharides formed by the cells in the liquid cultures, a process that continues after the cells are placed on the agarose plate. The formation of this \slime" facilitates motility and is necessary for the cells to move onto the non-pre-wetted surface. In (c) (as well as in (a) and (d)), the cells are observed to form small aggregations as a precursor to the nger of cells pushing the boundaries of wetted surface toward the direction of light. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 An example of using the Particle Tracker plugin for ImageJ. Cells that were identi ed in this frame have red circles around them. The light source imposes a force on the cells in the direction indicated on the picture. This same directional force is also present in the series of images in Fig. 2.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Example of trajectories produced using Particle Tracker. (a) Indi- cates cells which were stationary for the entire 2000 s. (c) Illustrates a pair of cells which move together for frames (a)-(e). (e) Highlights a cell which turns to move toward a neighbor. (g) Shows the direc- tion of light for the course of the experiment. Trajectories were only included if they were more than 100 frames (500 seconds) in length. Also, note that the di erent colors allow for optical di erentiation of cell trajectories and do not indicate anything about the cells or their motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 vii 2.4 An example of using Particle Tracker data to determine the direc- tion of the cells. In (a), we see a large mass of cells with a smaller region selected. This gure shows one image taken from a sequence of 301 with 2 seconds between each frame and spanning a total of 10 minutes. (b) Scatter plot showing the direction and how far each cell travels. (c) Histogram of the angular movement (in radians) between frames from the selected region in (a). The direction of light is to the upper left corner, or approximately 2.4 radians; however, note that the histogram in (c) does not indicate any preference of cells to move in that direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.5 (a){(c) Random motion model (Model R) and (d){(f) Random mo- tion with neighbor dependence model (Model B) on changing direc- tion for N = 250 cells at times t=0, 350, and 1000 time steps. The red circle around the cell in the upper left hand corner of each image represents the interaction distance D. For (a){(c), there is no rele- vant D and b = 0:5. For (d){(f), D = 5, a = 0:8 and b = 0:5. The particles are constrained to a circle with a radius of 100. . . . . . . . 26 2.6 (a) Partitioning the domain into 169 equal area sections. The num- bers represent partition number labels. (b) A histogram of the num- ber of cells in each partition over time for the random motion model (Model R). (c) A histogram of the number of cells in each parti- tion over time for the neighbor-dependent model (Model B) (d) A histogram of the number of cells in each partition over time for the local-interactions model (Model A) Each histogram includes N = 250 cells and counts the number of cells over t = 1000 time steps. For Model A and B, the parameters are D = 5, a = 0:8 and b = 0:5. For Model R, a = 0:8 and b = 0:5. Cells are clearly observed forming aggregations on the boundary in (c) and (d). . . . . . . . . . . . . . . 27 2.7 (a){(f) Model A; (g){(l) Model B. N = 250 Particles are constrained to a circle or diameter 200. Shown are simulation results at times t = 0; 50; 350. The interaction distance is D = 5 for (a){(c), (g){(i) and D = 10 for (d){(f),(j){(l). Identical initial conditions are used in all simulations as well as parameters a = 0:8 and b = 0:5. Note the di erent aggregation patterns and their dependence on parameter D. 29 2.8 (a), (b) Model A; (c), (d) Model B. Both models are simulated in R2 at t = 1000 for N = 250 particles. The interaction distance is D = 5 for (a), (c) and D = 10 for (b), (d). The initial conditions and parameters a, b are identical to those used in Fig. 2.7. More particles stay within the shown domain for simulations of Model A. . . . . . . 31 viii 2.9 (a), (b) Model A; (c), (d) Model B. Periodic boundary conditions; t = 1000; N = 250 particles. The interaction distance is D = 5 for (a), (c) andD = 10 for (b), (d). The initial conditions and parameters a, b are identical to those in Fig. 2.7. . . . . . . . . . . . . . . . . . . 33 2.10 Histogram illustrating how many cells have 0 to 10 neighbors within 5 or 10 pixels. Sub gures (a), (b) correspond to the simulation of Model A with D = 5 in Fig. 2.9(a). Sub gures (c), (d) correspond to the simulation of Model B with D = 5 in Fig. 2.9(c). In (a) and (c), we graph the number of neighbors that are within 5 pixels of each cell. In (b) and (d), we graph the number of neighbors that are within 10 pixels of each cell. . . . . . . . . . . . . . . . . . . . . . . 34 2.11 The same ve cells traced over 1000 time steps for (a), (b) Model A, and (c), (d) Model B. Periodic boundary conditions; t = 1000; N = 250 particles. The interaction distance is D = 5 for (a), (c) and D = 10 for (b), (d). The initial conditions and parameters a, b are identical to those in Fig. 2.7. Model A allows for more persistence in directed movement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.12 (a){(f) Six realizations of Model A with the same initial conditions and parameters a and b, using D = 10 and N = 250 cells at time t = 1000 time steps. While precise locations of aggregations vary, the aggregation sizes and numbers are similar. . . . . . . . . . . . . . . . 36 2.13 A study on the e ect of interaction distance D. We consider (a), (d) D = 10, (b), (e) D = 20, and (c), (f) D = 40 with the same initial conditions and same parameter a, b and N = 250 cells at time t=1000. For (a)-(c), the cells are constrained to a 200 pixel diameter circle. For (d)-(f), the cells are allowed to move anywhere in the R2 plane. As D increases, we observe fewer aggregations but with more particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1 (a){(c) Local interaction model and (d) Global forcing in direction of light. Model particles can exhibit (a) persistence, (b) stationary behavior, (c) movement toward neighboring cells, and (d) movement in the direction of a light source. . . . . . . . . . . . . . . . . . . . . 44 3.2 Local interaction Model A at t = 0; 50; 350 time steps with (a)-(c) D = 5, N = 250, (d)-(f) D = 10, N = 250, (g)-(i) D = 5, N = 500, and (j)-(l) D = 10, N = 500. . . . . . . . . . . . . . . . . . . . . . . 48 3.3 Global forcing Model H at t = 0; 500; 1000 time steps with (a)-(c) D = 5, N = 250, (d)-(f) D = 10, N = 250, (g)-(i) D = 5, N = 500, and (j)-(l) D = 10, N = 500. The direction of light is = =4. . . . 51 ix 3.4 Global forcing Model I with leader (red ) to normal (black ) cell ratio 1:1 at t = 0; 500; 1000 time steps with conditions (a)-(c) D = 5, N = 250, (d)-(f) D = 10, N = 250, (g)-(i) D = 5, N = 500, and (j)-(l) D = 10, N = 500. The direction of light is = =4. . . . . . 52 3.5 Long time simulations (t = 3000 and 9999) of each model for N = 250 particles with (a)-(f) D = 5 and (g)-(l) D = 10. The direction of light is = =4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.6 Simulations of global Model I at t = 0; 1000; and 3000 time steps for D = 5 with various ratios of leader (red ) to normal (black ) cells: (a)-(c) 1 leader : 3 normal, (d)-(f) 1 leader : 1 normal, (g)-(i) 3 leader : 1 normal and (j)-(l) All leader cells, i.e. global Model H. The direction of light is = =4. . . . . . . . . . . . . . . . . . . . 56 3.7 The activation model with (a){(c) K1 = 2, (d){(f) K1 = 4, and (g){ (i) K1 = 8 for excitation cuto K = 1:5. N = 250 particles are constrained to a circle of diameter 200. Shown are simulation results at times t = 50; 350; 1000. Identical initial conditions are used in all simulations as well as parameters D = 5, a = 0:8, c = 0:01 and b = 0:5. The direction of light is = 3 =4. . . . . . . . . . . . . . . 61 3.8 The activation model with at long-time t = 3000 for (a){(f) cuto K = 1:5 and (g){(l) K = 2:5. The detection distance is D = 5 in (a){ (c) and (g){(i) and D = 10 in (d){(f) and (j){(l). The rst column has K1 = 2, the second column has K1 = 4, and third column has K1 = 8. In each simulation, N = 250 particles are constrained to a circle of diameter 200. Identical initial conditions are used in all simulations as well as parameters a = 0:8, c = 0:01 and b = 0:5. The direction of light is = 3 =4. . . . . . . . . . . . . . . . . . . . . . 62 3.9 Scatter plot of the excitation levels at long-time t = 3000 for (a){(f) cuto K = 1:5 and (g){(l) K = 2:5. The y-axis is excitation value Si(t) and the x-axis is time t. In each gure, the excitations for all 250 cells are plotted. The detection distance is D = 5 in (a){(c) and (g){(i) and D = 10 in (d){(f) and (j){(l). The rst column has K1 = 2, the second column has K1 = 4, and third column has K1 = 8. In each simulation, N = 250 particles are constrained to a circle of diameter 200. Identical initial conditions are used in all simulations as well as parameters a = 0:8, c = 0:01 and b = 0:5. The direction of light is = 3 =4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 x 3.10 Examples of cell trajectories over t = 1000 time steps for the activa- tion model. The detection distance is D = 5 in (a){(c) and D = 10 in (d){(f). The rst column has K1 = 2, the second column has K1 = 4, and third column has K1 = 8. In each simulation, N = 250 particles are constrained to a circle of diameter 200. Identical initial conditions are used in all simulations as well as parameters K = 1:5, a = 0:8, c = 0:01 and b = 0:5. The direction of light is = 3 =4. . . . . . . . 65 4.1 A number line illustrating the cell position xi of particle i with de- tectable neighbors being the total number of particles which fall in the regions indicated by +i and i . In this particular image, the neighbor detection distance D = 6. . . . . . . . . . . . . . . . . . . . 70 4.2 A sketch of the notation used for describing the motions of particles. Not shown in this gure are the arrows showing particles in bin i 1 moving into population ri and showing particles in bin i + 1 moving into population li. a is the persistence probability and cri and c l i are the neighbor-weighted probabilities of a particle in bin i choosing to move to the right and left, respectively. . . . . . . . . . . . . . . . . . 73 4.3 Setup of particles for J+i;r[~r] . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4 Setup of particles for J+i;l[~r;~l] . . . . . . . . . . . . . . . . . . . . . . . 78 4.5 A depiction of the options of motion for particles out of a right-moving ri or left-moving li position, with associated probabilities for a system of particles which are either right-moving or left-moving and have the ability to become stationary with memory of a preferred direction. Not shown are the sources of right-moving and left-moving particles. The source of particles in ri is all particle types in bin i 1, and the source of particles in li is all particle types in bin i+1. The associated probabilities of these sources can be extracted from the diagram. The probability a denotes persistence, the probability b denotes becoming (or staying) stationary and the probabilities cri and c l i denote choosing a direction based on neighbor weights. . . . . . . . . . . . . . . . . . 83 4.6 A sample simulation of model (4.29){(4.32) with model parameters a = 0:5, b = 0:1 and D = 7. In (a), we show the evolution of the system in 100 bins over 1000 time steps. In (b), we show the state of the system at t = 1000. . . . . . . . . . . . . . . . . . . . . . . . . . 92 xi 4.7 A comparison of the stochastic particle system to the determinis- tic ODE for di erent scenarios: forming aggregations, the absence of aggregations, and merging aggregations. The parameters a, the persistence probability, D, the neighbor detection distance, and the initial data are identical in each case. The parameters a and D are varied as indicated above each plot to capture the di erent scenarios. For each of these images, the stopping probability b is 0.1 and the number of bins k is 100. On these 3-D images, the i-axis is for bin number, from 1 to 100, the t-axis is for time, from 1 to 1000, and the vertical z-axis is for particle number. . . . . . . . . . . . . . . . . . . 93 4.8 Comparison of the stochastic particle system and the deterministic ODE for di erent scenarios: aggregations, the absence of aggrega- tions, and merging aggregations at t = 1000. The parameters a, the persistence probability, D, the neighbor detection distance, and the initial data are the same in each case. The parameters a and D are varied as indicated above each plot to create the di erent scenarios. In all images, the stopping probability b is 0.1 and the number of slots k is 100. On these 2-D plots, the i-axis is for bin number, ranging from 1 to 100 and the vertical z-axis is for the number of particles, speci cally Mi(1000). . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.9 Four simulations of the stochastic model for the same parameter set, alongside the simulation produced by the deterministic ODE model. In all images, the persistence probability a is 0.5, the neighbor detec- tion distance D is 7, the stopping probability b is 0.1, and the number of bins k is 100. On the 3-D plots, the i-axis is for bin number, rang- ing from 1 to 100, the t-axis is for time step, from 1 to 1000, and the vertical z-axis is for total particle number, Mi(t). On the 2-D plots, the i-axis is for bin number and the vertical z-axis is for total particle number, Mi(1000). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.10 Exploration of parameter space for parameters a, the persistence probability, and b, the stopping probability, for uniform initial condi- tions (Mi(0) = 8 for all i). Note that there is a hard constraint that a+ b 1. For cases where a+ b = 1, particles do not have the ability to choose to move in a new direction; they can only stop or persist in their initially preferred direction. In all cases, the interaction dis- tance D is 10 and the number of bins k is 100. The i-axis is the bin number, ranging from 1 to 100, the t-axis is for time, from 1 to 1000, and the vertical z-axis is for total particle number, Mi(t). . . . . . . 99 xii 4.11 Exploration of parameter space for parameters a, the persistence probability, and b, the stopping probability, for uniform initial con- ditions (Mi(0) = 8 for all i). This gure has the same parameter setup and initial conditions as those in 4.10; however, the MATLAB integration tolerances are stricter in this Figure to illustrate the nu- merical instability of these simulations. Note that there is a hard constraint that a+ b 1. For cases where a+ b = 1, particles do not have the ability to choose to move in a new direction; they can only stop or persist in their initially preferred direction. In all cases, the interaction distance D is 10 and the number of bins k is 100. The i-axis is for bin number, ranging from 1 to 100, the t-axis is for time, from 1 to 1000, and the vertical z-axis is for particle number. . . . . 100 4.12 Initial conditions for Figures 4.13{4.22. The rst four plots, for right-moving Ri(0), right-stationary Rsi (0), left-moving Li(0) and left- stationary Lsi(0) are used to initiate the ODEs (and the stochastic systems in Fig. 4.7, 4.8 and 4.9. The last plot is the sum of these four plots, Mi(0) and is what is visible in subsequent images. Particle numbers are distributed with a Poisson distribution with mean 2 for each population-type in each bin i. The i-axis varies from 0 to 100. The total number of particles is 820. . . . . . . . . . . . . . . . . . . 101 4.13 Exploring the parameter space for parameters a, the persistence prob- ability, and b, the stopping probability. Note that a+b 1. For cases where a+b = 1, particles cannot change their direction; they can only stop or persist in their initially preferred direction. In all cases, the interaction distance D is 10 and the number of bins is 100. The i- axis is for the bin number, the t-axis is time, from 1 to 1000, and the vertical z-axis is for particle number. The initial conditions are given in Figure 4.12. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.14 The nal distribution at time t = 1000 of the simulations shown in Figure 4.13. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.15 Exploring the parameter space for parameters a, the persistence prob- ability, and D, the neighbor detection distance. In all images, the stopping probability b is 0.1 and the number of bins is 100. The i- axis is for the bin number, from 1 to 100, the t-axis is the time, from 1 to 1000, and the vertical z-axis is for the particle number. The initial conditions are given in Figure 4.12. . . . . . . . . . . . . . . . 106 4.16 The nal distribution at time t = 1000 of the simulations shown in Figure 4.15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 xiii 4.17 Exploring the parameter space for parameters a, the persistence prob- ability, and D, the neighbor detection distance. In all images, the stopping probability b is 0.1 and the number of bins is 100. The i- axis is for bin number, ranging from 1 to 100, the t-axis is the time, from 1 to 1000, and the vertical z-axis is for the particle number. The initial conditions are given in Figure 4.12. . . . . . . . . . . . . . . . 108 4.18 The nal distribution at time t = 1000 of the simulations shown in Figure 4.17. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.19 Exploring the parameter space for parameters b, the persistence prob- ability, and D, the neighbor detection distance. In all images, the stopping probability a is 0.3 and the number of bins is 100. The i- axis is for the bin number, from 1 to 100, the t-axis is for the time, from 1 to 1000, and the vertical z-axis is for the particle number. The initial conditions are given in Figure 4.12. . . . . . . . . . . . . . . . 111 4.20 The nal distribution at time t = 1000 of the simulations shown in Figure 4.19. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.21 Exploring the parameter space for parameters b, the persistence prob- ability, and D, the neighbor detection distance. In all images, the stopping probability a is 0.3 and the number of slots k is 100. The i-axis is for bin number, from 1 to 100, the t-axis is for time, from 1 to 1000, and the vertical z-axis is for the particle number. The initial conditions are given in Figure 4.12. . . . . . . . . . . . . . . . 113 4.22 The nal distribution at time t = 1000 of the simulations shown in Figure 4.21. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.1 Cartoon depiction of an e ector cell, or CTL, interacting with a can- cer cell via Fas/FasL binding to form complex X . . . . . . . . . . . 120 5.2 Percent lysis experimental data t by model at 4 and 12 hours for both B7-H1+ and B7-H1 for various e ector to target cell ratios. The experimental data is from [41]. . . . . . . . . . . . . . . . . . . . 128 5.3 Population dynamics for cancer cells, transfected with either B7-H1+ or a mock protein in the presence of e ector cells T and in the absence of e ector cells T . The concentration of the Fas-FasL complexes is also plotted but very small relative to T . For the initial conditions, we chose E=T = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 xiv 5.4 Latin hypercube sampling analysis of parameter k4, apoptosis rate for perforin interacting with cancer cells, where all other parameters are xed at values given in table 5.2. The value being plotted is the weighted sum of square errors for the plotted parameter values divided by the weighted sum of square errors for the parameter values in the table. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 xv List of Abbreviations 51Cr A radioactive isotope of Chromium with atomic weight 51 CTL cytotoxic T cell eq. equation FasL Fas ligand Fig. Figure ODE ordinary di erential equation PCC Pasteur Culture Collection RDME reaction-di usion master equation sp. species TLVM Time-lapse video microscopy xvi Chapter 1 Introduction Cells are very complex, microscopic, organic machines, composing all known life forms. They are able to convert multitudinous biochemical reactions into func- tional beings, from unicellular bacteria to multicellular animals. Despite the small size of cells, many open questions exist regarding cellular interactions and how such interactions produce observable macroscopic dynamics. Approaching these problems within a mathematical framework allows for a qualitative understanding of the un- derlying mechanisms. This dissertation is on mathematical models and simulations of intercellular interactions for two open research areas in biology and biomedical sciences: 1. Phototaxis and local interactions between cyanobacteria, and 2. Cancer-immune interactions and B7-H1, a surface protein found on some can- cer cells. The rst topic, phototaxis, is microbiological in nature and is the primary focus of this dissertation. It is discussed in Chapters 2{4. The second topic, which is related to cancer dynamics, is discussed in Chapter 5. In this introduction we highlight the biological and mathematical modeling background for both topics and provide an overview of the original work carried out in this dissertation. 1 1.1 Phototaxis 1.1.1 Biological background Unicellular microorganisms have evolved to live in variable and extreme envi- ronments. Some are capable of intercellular signaling and appear to utilize group dynamics to achieve desired actions, such as moving toward a food source [46] or to- wards light [11]. These group dynamics often result in emergent patterns which can be modeled and analyzed using mathematics. Plants, such as sun owers, can turn themselves to improve their ability to sense light and engage in photosynthesis, a phenomenon known as phototrophism. Phototaxis instead involves an entire organ- ism moving in response to light. Phototaxis allows for optimizing photosynthesis, habitat selection and improved reproductive ability [10, 37]. A variety of unicellular bacteria are phototactic, e.g. Chlamydomonas, Euglena, and Synechocystis [37, 9]. These di erent cells vary widely in size, from 1 micron to over 50 microns, and use a variety of di erent photoreceptors to sense light [37]. Some cells use information from two photoreceptors to determine the direction of a light source [50]; whereas some phototactic cyanobacteria use one photoreceptor to sense light at di erent time points, using a di erencing method to identify the direction of light [39]. The process of identifying the direction of light is not known for all species, e.g. Synechocystis sp. Phototactic cells take many physical forms; some phototactic cells are agel- lar, some are ciliates, and some have pili. Di erent cell motility structures result in di erent global movement patterns; motility produced by these various cellular structures is studied by biologists and mathematicians alike [23, 60, 62, 76, 80]. 2 We are primarily concerned with the common freshwater cyanobacterium Syne- chocystis sp. Synechocystis species PCC6803 (hereafter Synechocystis sp.) displays the ability to move toward light, forming nger-like projections in the direction of a light source. Synechocystis sp. is a model organism for studying phototaxis in a lab- oratory setting. Extensive genetic and microscopic analyses have been carried out to characterize the molecular basis for motility [9]. It has been demonstrated that this surface dependent motility requires Type IV pili, unusual photoreceptors [44] and a host of other less well-characterized proteins [9]. The pili are long hairlike ap- pendages of di erent diameters. Some pili are used as the motion apparatus, while the role of other pili is lesser known, though it is likely that they play some role in signaling between bacteria. The movement and pattern formations have been experimentally shown to depend on the density of cells, wavelength of light, and the properties of the surface [11]. Phototaxis in Synechocystis sp. is a complex and well-regulated process. Several mutants in phototaxis have been identi ed [9] and yet a complete understanding of the signals and intercellular interactions that con- trol the emergent behaviors of phototactic cells is still lacking. When Synechocystis sp. is exposed to light, some cells begin to move, although not immediately and not necessarily in the direction of light. They instead form small aggregations of cells and eventually, with a time delay, some cells may move toward light, forming the characteristic nger-like projections or swarms of cells [15]. Yet the motion of individual cells is not as directed toward the light source as is the observed group behavior. Individual cells instead display a quasi-random motion, that is, they move in seemingly random directions. 3 1.1.2 Mathematical models relevant to phototaxis An extensive amount of mathematical research has been conducted in related elds. Speci cally, we would like to mention some works on animal ocking and chemotaxis. The Couzin-Vicsek model of ocking (and its many extensions) allows animals or individual agents to be repelled by near neighbors, align with the average directional heading of not-so-near neighbors, and be attracted to far neighbors [17, 79]. This model has lent itself to many applications as well as thorough mathematical analysis, for example see [22]. Cucker and Smale o er a dynamical system which models the ocking of opinions in human networks [18, 19]; this model has also been subjected to signi cant analysis, for example see [35, 36]. Similar models of ocks and schools have been developed for a variety of self-propelling agents such as birds and sh, e.g. [3, 43, 57, 58, 66]. Chemotaxis, i.e., motion towards a chemical attractant, is also a eld that has been extensively studied in recent decades, starting from the celebrated work of Patlak, Keller and Segel [47, 67]. For completeness, we refer the interested reader to the following papers and to the references therein [1, 40, 42, 65, 78]. The motion of chemotactic cells is typically subjected to di erent rules of motion than those that are obeyed by Synechocystis sp. Consequently, most of the mathematical modeling and analysis on the topic is irrelevant in the present context. Furthermore, phototaxis has not been extensively subjected to mathematical modeling. Only a few models of phototaxis have been developed, for example see [14, 59]; however, these models do not consider the intercellular group dynamics. A 4 recent agent-based model of phototaxis considers cell interactions by transmission of light by individual cells [30]. Another recent ODE and statistical model of pho- totaxis examines rotational properties of an algal colony of bi agellar V. Carteri cells [26]. In a series of papers, my collaborators have developed several families of mathematical models for describing phototaxis [13, 15, 54, 55, 56]. In all these models, the primary focus was on the phase of the initiation of the motion towards a light source (including the associated time delay) and the resulting overall mi- gration of the colony of cells towards the light source (including the modeling of the nger formation). The emphasis of these analyses was on the role of the group dynamics, as opposed to what can be associated with the behavior of the individual cell. Missing from these analyses was the description of what happens in regions of low to medium cell density. In this our work, we seek to develop, present, and study mathematical models for such regions and that account for experimentally observed quasi-random motion of cyanobacteria. 1.2 B7-H1 and tumor-immunology 1.2.1 Biological background B7-H1 is a surface protein which has been found on carcinomas of the lung, ovary and colon and in melanomas but not on most normal tissues [25]. An un- derstanding of a blockade of B7-H1 has been theorized to be applicable not only to cancer, but also to viral infections and in ammatory bowel disease [16]. In the case of tumor immunity, experiments have shown that when the B7-H1 surface pro- 5 tein, also referred to as PD-L1 and CD274, is present on a tumor, cytotoxic T cells (CTLs) are less e ective at inducing apoptosis in the cancer cells [25, 41]. It is believed that B7-H1 forms a blockade against CTLs by interacting with PD-1 on the surface of CTLs [4]. Phase I clinical trials with an antibody which blocks B7- H1/PD-1 interactions yielded complete remission in one patient with non-Hodgkin?s lymphoma [8]. A better understanding of the mechanism and dynamics may allow medical researchers to develop a cancer treatment schedule speci cally targeting this molecular shield, allowing the immune system to more e ectively repress a tumor. Cytotoxic T cells interact with cancer cells as to induce apoptosis via two mechanisms: Fas/FasL binding and perforin [82, 2]. In the case of Fas/FasL binding, the Fas ligand (FasL), which is present on some CTLs, forms a complex with Fas, a protein on the surface of a tumor cell. Through this complex, CTLs are able to send apoptotic signals to the cancer cells. In the case of perforin, the CTLs produce perforin which is thought to perforate the tumor cell surface. This allows CTLs to inject cancer cells with granzymes which induce apoptosis. This overview of the apoptosis mechanism is clearly oversimpli ed. For a more thorough biological review of these mechanisms we refer to [70]. The Fas-FasL mechanism is theorized to be a ected by the presence of B7-H1 [25], while perforin-mediated lysis, the primary means of lysis during the rst four hours of culture, is expected to be less a ected by the presence of B7-H1 [41]. Cytotoxicity can be measured experimentally by performing a Chromium release assay which measures percent lysis of target, or cancer, cells by e ector cells, or CTLs. 6 1.2.2 Models relevant to B7-H1 and cancer immunology Many tumor-immune system interaction models have been developed to date. None of the existing mathematical models deal directly with B7-H1. A non-comprehensive list of references include the following. In 1994, Kuznetsov et al. developed a basic interaction model which utilizes experimental data to determine rate parameters. Kirschner and Panetta o er a tumor-immune interaction model which incorporates IL-2 [48]. DePillis et al. o er an ordinary di erential equation model which utilizes \Hill"-like terms and ts percent lysis data [24]. Thorn and Henney produced a rel- atively simple model using Michaelis-Menten enzyme-substrate kinetics which also utilizes percent lysis data [77]. Wilson and Levy use an ODE model, t to experi- mental data, to analyze the interplay of a vaccine and immunotherapy [83]. Models using kinetic theory have even been developed [6, 20, 21]. There are a few very mech- anistic models of CTL-induced apoptosis. Lai o ers a model of the Fas-FasL trimer complex formation biology [52]. Another model considers the e ect of Fas/FasL binding on a population scale [81]. A detailed model for the internal workings of granzyme-induced apoptosis has also been produced [34]. However, none of these models consider both mechanisms of CTL induced apoptosis of a target cell. The only mathematical modeling work to date on tumor-immune interactions speci cally addressing B7-H1 is our paper [33], presented in Chapter 5. 7 1.3 Making the connection between phototaxis and cancer Both phototaxis and cancer-immune interactions are characterized by inter- cellular interactions and communications which are not thoroughly understood, and both are amenable to mathematical modeling e orts. As research in these elds progresses, a few other links between these elds may emerge. One of the long-term applications of studying phototaxis is in light-directed medicine. Understanding how certain cells move under the in uence of light is critical to future technological advancements that may allow control of this movement. We might then be able to control cells, with drugs attached or genetic material encoded, to target speci c cancerous tissues using light. Certain related e orts are currently being undertaken using magnetic elds and drug-coated magnetic foreign bodies [61, 74], as well as light-activated drugs [7, 73]. The study of phototaxis may enhance our knowledge of such light-activatable molecules, which could be utilized to ensure drugs are only activated in the target area, so as to not endanger healthy tissue [38]. Furthermore, chemotaxis, the phenomenon of cells moving toward a chemoattractant, which shares some common principles with phototaxis, plays a fundamental role in cancer growth [68, 53, 69, 71, 75], and in immunology in general, e.g. [28, 49, 64, 84]. 1.4 Overview of dissertation In this dissertation, we develop mathematical models and simulations for two biological topics in intercellular interactions: phototaxis and B7-H1 cancer-immune interactions. 8 In our models of phototaxis, we seek to develop a better understanding of the local interactions and global forcing e ects of directed light on the phototactic cyanobacterium Synechocystis sp. We begin in Chapter 2, by analyzing experi- mental data which are used to develop two agent-based stochastic models of local interactions between particles. Numerous simulations are conducted to shed light on the role and impact of the model parameters. In Chapter 3, we extend our agent- based model from Chapter 2 to include global forcing due to a light source. This model accounts for the activation of motion of cells triggered by light. Such external forcing leads to aggregations forming and moving toward a light source with their new awareness of the direction. In Chapter 4, we reduce our original model of local interactions to a one-dimensional discrete stochastic particle system. This is used to derive a reaction-di usion master equation (RDME) for the probabilistic evolution of the system, from which we can derive a system of ordinary di erential equations (ODEs) which deterministically describe the motion due to local interactions on a line. This system of ODEs provides us with a direct and e cient access to study the parameter space. The simulations of our mathematical models qualitatively replicate many ex- perimentally observed behaviors of Synechocystis sp. Our modeling e orts support the novel idea that not all cyanobacteria are capable of sensing the direction of a light source, i.e., that di erent traits of this rather primitive organism may exist. In our model of B7-H1 and cancer-immune interactions, we study the e ect of the presence of surface protein B7-H1 on cancer-immune interactions. In Chapter 5, we develop a model for cancer-immune interactions which depends on the two mech- 9 anisms of apoptosis: perforin-mediated and Fas/Fas-L-mediated. We isolate a set of parameters which depend on the presence of surface protein B7-H1. We t our model to in vitro percent lysis data and support the predicted e ect of B7-H1 on the two mechanisms of apoptosis. Finally, in Chapter 6, we discuss general conclusions of the research contained in this dissertation. 10 Chapter 2 Modeling local interactions during motion of cyanobacteria 2.1 Introduction In this chapter, we mathematically model the local interactions of phototactic cyanobacterium Synechocystis sp. and observed quasi-random motion in order to address a series of questions: 1. If motile cells are not moving exclusively toward the light, are they moving in random directions? 2. Is there a characteristic distance, within which cells are able to sense the presence and behavior of neighbors? 3. If such a distance exists, what are the mathematical e ects and biological implications of varying this distance? 4. Can we develop a model which explains the presence of the observed aggrega- tions of cells? Our mathematical models follow the time-discrete dynamics of a nite set of particles that are interacting in a two-dimensional domain according to rules that involve certain random terms. The rules for the local interactions between particles are based on our experimental observations given by the analysis of time- lapse movies of the bacteria under a plethora of controlled conditions. 11 The structure of this chapter is as follows. In Section 2.2 we present the biological background, describe the experimental setup, demonstrate some of the experimental results, and formulate a list of observations that lead us in devel- oping the mathematical models. The mathematical models are then presented in Section 2.3. Two mathematical models are described. Our rst model is a discrete- time model of local interactions that allows cells to keep moving in their previous direction, stop moving, or move in the direction of one of their neighbors. The second model increases the randomness in the motion. This time, cells can elect to change their direction of motion, but only if there are other cells nearby. Simulations of these models are then presented in Section 2.4. These simulations show that the models, in particular the local-interactions model, provide results that replicate the experimental observations. Concluding results are given in Section 2.5. 2.2 Biological background In our study, we consider a phototactic microorganism Synechocystis sp. This unicellular cyanobacterium has been studied extensively by time-lapse video mi- croscopy (TLVM). The experimental observations that serve as the basis to the mathematical models that follow were generated in the Bhaya lab by our collabo- rators, Devaki Bhaya and Susanne Wisen, at the Carnegie Institute for Science at Stanford University. Time-lapse movies of moving bacteria were acquired under an optical microscope and then analyzed as described below. In Figure 2.1, we show several examples of some of the experimentally ob- 12 (a) (b) (c) (d) Figure 2.1: Experimental results illustrating (a) beginnings of nger-like projec- tion formation, (b) surface dependence, and (c/d) aggregation and nger formation. The bacteria within a matrix of extracellular polysaccharides appear as white dots or circles. The bacteria are approximately 1 m in diameter. Each image contains wildtype cells exposed to a light source radiating from the upper left corner of each image. Inset (a) shows the front of a spot of cells 48 hours after plating. Inset (c) shows the exact same spot frame as (a), 24 hours later. Images (b) and (d) depict di erent sections of the same spot at di erent times. In (b), the density of cells is very low and surface dependence is observable|the cells are surrounded by polysaccharides formed by the cells in the liquid cultures, a process that continues after the cells are placed on the agarose plate. The formation of this \slime" facili- tates motility and is necessary for the cells to move onto the non-pre-wetted surface. In (c) (as well as in (a) and (d)), the cells are observed to form small aggregations as a precursor to the nger of cells pushing the boundaries of wetted surface toward the direction of light. 13 served characteristics of motion. Each image contains wildtype Synechocystis sp., exposed to a directional light source coming from the top left of the domain. In Fig. 2.1(a), the beginnings of nger formation are observed. In Fig. 2.1(c), the same part of the spot is observed 24 hours after (a) with signi cantly more prominent ngers. The cells observed in Fig. 2.1(b) and 2.1(d) are taken from the same ex- periment at di erent time points and locations on the plate and highlight the wide range of observable cell behavior. Fig. 2.1(b) illustrates the e ect of the surface on movement { cells preferentially move on a pre-wetted surface which contains extra- cellular polysaccharides produced by the cells. In (d), cells form small aggregations and move toward the light-source. Such aggregations of cells can be seen in all four images. We also observe the tendency of cells to align along the boundary between the pre-wetted and agarose-only surface in Fig. 2.1(a), (c), and (d). In each gure, we note the tendency of some bacteria to move in pairs, although it is important to note that in some cases, these pairs of cells are cells which have divided but have not separated. Not illustrated in Fig. 2.1, the wavelength of light and density of cells also impact nger formation and general cellular movement. Additionally, directional cell movement is typically delayed; i.e. after cells are plated, they do not appear to start moving in the direction of light until a few hours have passed. Furthermore, cells have been observed to sense the presence of and move toward other cells over a large distance, as is shown in [55]. This seems to imply the existence of at least one communication mechanism between cells. To analyze the experimental data, we used ImageJ, a public domain image 14 processing software produced by the National Institutes of Health, and a plug-in called Particle Tracker developed by Sbalzarini and Koumoutsakos [72]. In Fig. 2.2, we demonstrate the cell-position-identifying capabilities of Particle Tracker. We observe that most of the cells in the image are correctly identi ed. @ @I @ @I @ @I to light Force due Figure 2.2: An example of using the Particle Tracker plugin for ImageJ. Cells that were identi ed in this frame have red circles around them. The light source imposes a force on the cells in the direction indicated on the picture. This same directional force is also present in the series of images in Fig. 2.3. After the cells have been identi ed in a series of images, we further use the plug-in to track cells as they move. In Fig. 2.3, the di erent sub gures are shown at consecutive times, taken 50 frames, or 100 seconds, apart. These cells, observed under an optical microscope, are exposed to a directional light source, the force pointing toward the upper left corner of the each image, as illustrated in Fig. 2.2 and Fig. 2.3(g). Particle Tracker is used to identify trajectories over 453 images, taken at 5 second intervals. In Figure 2.3, we show the 44 trajectories which are at least 100 frames in length. In Figure 2.3, a few patterns of motion can be identi ed using the trajectory analysis. In Fig. 2.3(a), a couple of cells are identi ed as stationary. The small dot-like trajectories of these cells illustrate that these cells do not move signi cantly 15 (a) t = 0 s (b) t = 250 s (c) t = 500 s stationary C CO n move as pair Z Z} (d) t = 750 s (e) t = 1000 s (f) t = 1250 s move toward neighborHj (g) t = 1500 s (h) t = 1750 s (i) t = 2000 s @@I @@I@@I to light Force due Figure 2.3: Example of trajectories produced using Particle Tracker. (a) Indicates cells which were stationary for the entire 2000 s. (c) Illustrates a pair of cells which move together for frames (a)-(e). (e) Highlights a cell which turns to move toward a neighbor. (g) Shows the direction of light for the course of the experiment. Trajectories were only included if they were more than 100 frames (500 seconds) in length. Also, note that the di erent colors allow for optical di erentiation of cell trajectories and do not indicate anything about the cells or their motion. 16 from their original positions. Cells moving as a pair are marked in Fig. 2.3(c). In Fig. 2.3(e), we see a cell that abruptly changes its direction to move toward a neighboring cell. Additionally, in comparing the series of gures, it is possible to identify at least one cell moving in the direction of light, although this is clearly not the case for the majority of cells. It is clear from Fig. 2.3 that the particles do not obey a simple rule of motion. They do not all move towards the light source. Is there any bias in the direction of their motion? To study this question we use numerical data generated by Particle Tracker to compute a histogram illustrating the distribution of angular direction of cellular movement. The data set we considered this time is shown in Fig. 2.4(a). We cropped this picture to the small rectangle shown in the northeast corner of the picture and used Particle Tracker to detect and track particles. The light source in this picture is to the northwest, or 34 radians (about 2.4), and directly above and to the left of the small selection, in the region with a visibly higher cellular density, is a nger moving toward the light. A scatter plot of the distances and di- rections travelled by the particles from the selected rectangle is shown in Fig. 2.4(b). Fig. 2.4(c) is a histogram for angular directional movement of cells between time frames. One might expect that the cells in the selected area would also be traveling in the direction of the nger, however, as can be seen in the histogram in Fig. 2.4(c) this is not the case. The cellular movement is almost uniformly distributed with no apparent bias toward light and possibly a slight bias toward moving to the right during this 10 minute interval. The focus of this work is on regions of low to medium density, where we observe 17 (a) (b) (c) Figure 2.4: An example of using Particle Tracker data to determine the direction of the cells. In (a), we see a large mass of cells with a smaller region selected. This gure shows one image taken from a sequence of 301 with 2 seconds between each frame and spanning a total of 10 minutes. (b) Scatter plot showing the direction and how far each cell travels. (c) Histogram of the angular movement (in radians) between frames from the selected region in (a). The direction of light is to the upper left corner, or approximately 2.4 radians; however, note that the histogram in (c) does not indicate any preference of cells to move in that direction. 18 a quasi-random motion. Cells do not necessarily move only towards light, but also in a di erent, seemingly random motion as demonstrated by the particle trajectories in Fig. 2.3. Even at the base of a directed nger as seen in Figure 2.4, the overall movement of cells is not directed toward the light source. This quasi-random motion has also been observed within ngers on short time scales. The basic characteristics of the observed dynamics of the cells in the low to medium density area can be described follows: 1. At any time, cells that are moving may stop moving. Conversely, at any time, stationary cells may start moving. This is not clear from the tracking example in Fig. 2.3, but has been observed in movies. 2. There is a clear persistence in the direction of the movement. Cells that move in a given direction may continue moving in the same direction. This is apparent in Fig. 2.3 with the relatively straight tracks that the cells follow. 3. The duration of the time in which a cell persists moving in a given direction is not xed. It can vary over a large range of values. This is also evident from Fig. 2.3 as the cells travel di erent distances before turning. 4. The direction of the motion of a cell may change at any time. While this is illustrated in Fig. 2.3, it is more obvious in the TLVM movies. 5. When the direction of the motion changes, there is some tendency of moving in the direction of a neighboring cell. One example of this was highlighted in Fig. 2.3. 19 While analysis of time-lapse movies and snapshots of Synechocystis sp. is useful to characterize general movement of cells, it is still quite di cult to understand the fundamental interactions that are taking place between cells that allows for the observed movements and patterns. 2.3 A mathematical model for local interactions of Synechocystis sp. In this section we present mathematical models of the quasi-random motion of cells. Two models are discussed: a local-interactions model and a neighbor- dependent random motion model. In exploring local interactions between cells, it is important to note that simple random motion does not elucidate the patterns we observe experimentally. This is evident, e.g., by the nature of the particle trajectories shown in Fig. 2.3. 2.3.1 Model assumptions We begin by stating our formative assumptions. 1. Neighbor detection. We assume that cells not only can detect the presence of neighbors within a given radius but also can detect the relative location of neighbors within that radius. The radius of a neighbor detection will be considered to be a model parameter. This assumption is motivated by a few biological observations. First, the pili which grow on the cell are known to get longer with time. Cells have been observed to make connections through these pili. These connections are not easily observed with an optical micro- 20 scope and are, hence, not shown in our images illustrating typical cell behavior in Fig 2.1. The length of the pili could be directly related to the neighbor de- tection distance. Second, a signaling molecule, which quickly degrades, could have a similar e ect. A larger neighbor detection distance could be observed over longer time scales or for di ering agarose densities. Lower densities may allow for faster di usion of a signaling molecule. Such a signaling molecule has not yet been experimentally identi ed for this bacterium, but it is reasonable to assume one may exist. Overall, not that much is known about either mech- anism biologically, and therefore we assume the existence of a nite neighbor detection distance without making any further assumptions about the under- lying biological mechanism, leading to our next assumption. 2. Communication mechanism unknown. As mentioned, we do not make any assumptions about the speci cs of the communication mechanisms which allows cells to locate neighbors; though various mechanisms are feasible, in- cluding chemical signaling, local force sensing, and physical interaction. While these mechanisms are of interest, they will not play any role in our model, at least at this stage. 3. Cell density. We also assume that cell density is relatively low and that spatial hindrance e ects are negligible so that we can use a particle model. 4. Movement. Regarding the motion of cells, we assume cells may follow one of the following rules: 21 (i) A cell may continue moving in the direction it was previously moving, (ii) A cell may stop moving, or (iii) A cell may orient itself to move in a di erent direction. For our local interaction model, we will assume that the new direction is set toward a neighboring cell. For the two random models discussed above, any direction is equally viable. 5. Memory. One of the key assumptions embedded in the model is that cells have very long memory (in fact, the way the rules are formulated, the memory lasts forever). This means, for example, that a cell that became stationary still remembers the direction of the motion that brought it to its state, and hence motion in the same direction can resume at a later point in time. In such a way, particles that aggregate and stick together can nd themselves separated at a later time. While this has not yet been shown for Synechocystis sp., the Myxococcus bacterium has been observed to leave a polarized trail of slime, detectable by the cell, which supports the idea of a cell having a long-lasting memory [85]. 6. No directional persistence with change of direction. We also observe that, unlike agellar chemotactic cells, Synechocystis sp. have no visually apparent directionality. That is, the cells do not have an obvious head or tail which indicates the preferred directionality of the cell and whether the cell is moving forward or turning. We model this by assuming that if a cell turns, the turning angle is not biased based on the previous direction. Other 22 mathematical models have included a bias in which agents, when turning, tend to choose new directions that are similar to the current forward-moving direction, e.g. [63]. That is, we assume our cells have no inherent directionality as we do not have experimental evidence to justify imposing such a bias, especially in the phase of cells at the base of a nger where quasi-random motion is prevalent. Note that this type of bias is di erent from considering global forcing on the cells due to a directed light source. 2.3.2 Model A: A local-interactions model Given the three rules of motion, we formulate our model as a discrete-time dy- namical particle system. Our approach resembles the method used to study persis- tence in chemotaxis in [63]. The state of the system at time t is fxi(t); vi(t); i(t)gNi=1 where N is the number of cells, vi(t) 2 f0; 1g is the velocity of cell i, i(t) 2 [0; 2 ) is the angular direction of cell i, and xi(t) 2 R2 is the position of cell i calculated using vi(t), i(t) and xi(t 1). Note that the velocity of each cell can be either 1 or 0. In this way, we are assuming that a cell can either move at constant velocity or stop, as discussed in the assumptions above. In this model, we allow cells to move in the direction of a neighboring cell. We de ne a neighboring cell to be any cell j which is within an interaction distance D from the cell i of interest. In this stochastic model, we only allow a cell to change 23 direction with probability 1 a. We formulate this as: i(t+ 1) = 8 >>>< >>>: i(t); with probability a, j; with probability 1 ani for j s.t. xj 2 BD(xi), (2.1) where BD(xi) is a ball of radius D centered at xi, ni is the number of cells with position xj contained in BD(xi), i.e. jjxi(t) xj(t)jj2 < D, and j is the angle of the vector pointing from particle i to particle j. In this way, when a cell changes direction, with probability 1 a, the direction presented by each neighboring cell j has an an equal probability, (1 a)=ni, of being selected as the new directional heading for cell i. Note that if ni = 0, a cell does not have any neighbors. In this case the variable i remains unchanged, i.e., i(t+ 1) = i(t). The velocity in our model (which is either 0 or 1) is determined independently of the choice of direction. It remains unchanged with probability b. It switches to the other state with probability 1 b. This is formulated as vi(t+ 1) = 8 >>>< >>>: vi(t); with probability b, 1 vi(t); with probability 1 b. (2.2) Finally, after updating the velocity and angle associated with each cell, all cell positions are updated according to xi(t+ 1) = xi(t) + vi(t+ 1) 2 6 6 4 cos( i(t+ 1)) sin( i(t+ 1)) 3 7 7 5 : (2.3) 24 2.3.3 Model B: Random neighbor-dependent movement model In the random neighbor-dependent movement model, a cell is only allowed to change direction if there is a neighbor within an interaction distance D. In this case, the direction is chosen randomly. We formulate this as: i(t+ 1) = 8 >>>< >>>: i(t); with probability a, (t); with probability 1 a if there exists xj 2 BD(xi). (2.4) Here, for any time t, (t) is a uniform random variable distributed in [0; 2 ). The position and velocity are calculated in the same way as in Section 2.3.2. 2.4 Simulations To justify the complexity of the Models A and B, let us rst consider a simpler model. Consider a discrete-time model of random motion where cells are allowed to change direction to a uniformly chosen random variable on [0, 2 ) with probability 0.2 for each time step; we refer to this as Model R. We compare this random motion Model R with a the neighbor-dependent random model, Model B, where, for each time step with probability 0.2, cells are allowed to change direction to an angle chosen uniformly from [0, 2 ), but only if there is a neighboring cell within 5 pixels distance. In Figure 2.5, we compare the random model, shown in Fig. 2.5(a){(c), to Model B, where a directional change in movement depends on the presence of a neighbor within D = 5 pixels. Clearly, Model B, which is a rather simple neighbor- 25 dependent model, yields results that better match the experimental data, compared with the random model. This is manifested in the cells in Fig. 2.5(d){(f) that tend to align along the boundary (similar to what is experimentally observed). (a) (b) (c) (d) (e) (f) Figure 2.5: (a){(c) Random motion model (Model R) and (d){(f) Random motion with neighbor dependence model (Model B) on changing direction for N = 250 cells at times t=0, 350, and 1000 time steps. The red circle around the cell in the upper left hand corner of each image represents the interaction distance D. For (a){(c), there is no relevant D and b = 0:5. For (d){(f), D = 5, a = 0:8 and b = 0:5. The particles are constrained to a circle with a radius of 100. We quantify this behavior in Figure 2.6. To do this, we partitioned the 100 pixel radius circle into 169 equal area sections. The sections are labeled in Fig. 2.6(a). For each model, we then identi ed the numbered partition for each cell in space and times from start to 1000 time steps. This data was then compiled into a histogram counting how many cells fall into each section for each model. In this way, each cell is counted 1001 times. In Fig. 2.6(b), we show a histogram of the cell locations for the random motion Model R. We observe that the cells are relatively uniformly 26 distributed in space-time. In contrast, in Fig. 2.6(c), we show a histogram of the cell locations over each time step for the neighbor-dependent random motion, Model B. In this case, we observe that cells are gathering on the boundary, partitions 122 through 169, with approximately twice the probability of them falling in the interior of the circle, partitions 1 through 121. This ratio is likely parameter dependent. Qualitatively, this agrees with the experimental observations. (a) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108109 110111112113114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130131132133134135136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154155 156157158159160161 162 163 164 165 166 167 168 169 (b) 0 50 100 150 0 1000 2000 3000 4000 (c) 0 50 100 150 0 1000 2000 3000 4000 (d) 0 50 100 150 0 1000 2000 3000 4000 Figure 2.6: (a) Partitioning the domain into 169 equal area sections. The numbers represent partition number labels. (b) A histogram of the number of cells in each partition over time for the random motion model (Model R). (c) A histogram of the number of cells in each partition over time for the neighbor-dependent model (Model B) (d) A histogram of the number of cells in each partition over time for the local-interactions model (Model A) Each histogram includes N = 250 cells and counts the number of cells over t = 1000 time steps. For Model A and B, the parameters are D = 5, a = 0:8 and b = 0:5. For Model R, a = 0:8 and b = 0:5. Cells are clearly observed forming aggregations on the boundary in (c) and (d). 27 We now consider a few simulations comparing Model A, the local-interactions model, from Section 2.3.2, to Model B, the neighbor-dependent random movement model, from Section 2.3.3. In Fig. 2.7, we see a comparison of Models A and B in which cells are con- strained to a circle with a diameter of 200. We compare the e ect of varying inter- action distance and consider the cases of D = 5 and D = 10. The length scale of interaction distance D is shown in the upper left corner of each picture around a stationary cell. This cell, contained in a circle of radius D, is not included in any of the simulations. Each simulation is run with N = 250 cells and can be seen at t = 0; 50; and 350 increments. The initial conditions, seen in 2.7(a), 2.7(d), 2.7(g), and 2.7(j), are identical in all simulations. We observe that for the local interaction model in 2.7(c) and the random neighbor-dependent motion model in 2.7(i), the same cells tend to stay around the border, similarly to what is seen in experiments. We recall that there is no mechanism that is embedded into the model that forces cells to stay on the boundary of the domain once the boundary is reached. Indeed, occasionally, individual particles return back into the domain. The cells appear to be slightly more uniformly distributed in 2.7(i) than in 2.7(c). In 2.7(f), Model A with D = 10 yields cells aggregating in small clumps of 5 to 15 cells. The simu- lated cells in Fig. 2.7(l) appear to be distributed uniformly with little cohesion to the boundary. In fact, the motion observed in the simulation depicted in 2.7(j){(l) is very similar to the random motion model from Figure 2.5 because the chosen surface density and interaction distance are such that, on average, each cell has 1.5 neighbors, e ectively removing the neighbor-dependence from this model. 28 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 2.7: (a){(f) Model A; (g){(l) Model B. N = 250 Particles are constrained to a circle or diameter 200. Shown are simulation results at times t = 0; 50; 350. The interaction distance is D = 5 for (a){(c), (g){(i) and D = 10 for (d){(f),(j){ (l). Identical initial conditions are used in all simulations as well as parameters a = 0:8 and b = 0:5. Note the di erent aggregation patterns and their dependence on parameter D. 29 In Figure 2.6, we o er a quantitative comparison of Models A and B for D = 5. We partition the domain into 169 equal-area sections and count the total number of particles that are present in each cell over 1000 time steps. We observe that for Model A, in Fig. 2.6(d), the boundary, partitions 122 through 169, has more cells on it with larger variability than for Model B (shown in Fig 2.6(c)). We suspect that this is due to the formation of aggregations on the boundary. As previously stated, cells have been experimentally observed to collect along the boundary. The boundary conditions in Fig. 2.7 simulate a setup that is similar to the experimental setup focusing on short-time dynamics of cells. In experiments lasting relatively short periods of time (on the order of hours), cells are con ned to a given area due to the properties of the underlying surface. However, in experiments lasting longer periods of time (on the order of days), cells overcome the underlying surface and, in the absence of a directional global light source, spread out on the surface of the agarose. In contrast, we consider large-time simulations of cells that are free to move anywhere in R2. Such simulations are shown in Fig. 2.8. The initial conditions are identical to the initial circular drop of cells seen in Fig. 2.8(a). The local- interactions model, Model A, is shown in Fig. 2.8(a) with D = 5 and Fig. 2.8(b) with D = 10. The random neighbor-dependent model, Model B, is shown in Fig. 2.8(c) with D = 5 and Fig. 2.8(d) with D = 10. All results are shown at t = 1000. Clearly, in Figs. 2.8(a),(c),(d), the majority of cells have dispersed. Fig. 2.8(b) stands di erent. In this case, of Model A, with the longer interaction distance, few aggregations still remain after these rather long-time simulations, in spite of the 30 particles not being con ned to the domain that is shown in the gure. D = 5 D = 10 (a) (b) M ode l A (c) (d) M ode l B Figure 2.8: (a), (b) Model A; (c), (d) Model B. Both models are simulated in R2 at t = 1000 for N = 250 particles. The interaction distance is D = 5 for (a), (c) and D = 10 for (b), (d). The initial conditions and parameters a, b are identical to those used in Fig. 2.7. More particles stay within the shown domain for simulations of Model A. To complement the previous cases, we also consider periodic boundary condi- tions. In Fig. 2.9, we show results obtained from simulations of cells on a 240 240 re- gion with periodic boundary conditions. The same initial conditions as in Fig. 2.7(a) are used. The results are shown at time t = 3000. Model A is shown in Fig. 2.9(a) with D = 5 and in Fig. 2.9(b) with D = 10. Model B is shown in Fig. 2.9(c) with D = 5 and in Fig. 2.9(d) with D = 10. In these snapshots, aggregations are visible in Fig. 2.9(b) for D = 10 in Model A. There appears to be slightly more white space in Fig. 2.9(a) than in Fig. 2.9(c) and Fig. 2.9(d). In order to quantitatively com- 31 pare this clumping e ect seen in Fig. 2.9(a) and Fig. 2.9(c), we show a histogram of how many neighbors are found within a distance of 5 or 10 pixels of each cell in Fig. 2.10. For Model A with D = 5, we observe in Fig. 2.10(a) that within 5 pixels of each cell, there are more cells with two or more neighbors than in Model B shown in Fig. 2.10(c). In comparing the same histograms, there are also slightly more cells with no neighbors for Model B than for Model A. Within 10 pixels, the phenomenon is slightly reversed; that is, within 10 pixels, in Fig. 2.10(b) there are more cells with no neighbors for Model A than for Model B, in Fig. 2.10(d). This trend reversal implies that the aggregates seen in simulations of Model A may be tighter, as there are more isolated cells, meaning the other cells must be packed into a slightly smaller space. In Fig. 2.11, the trajectories of ve cells, with identical initial conditions and the same model and interaction distance conditions as in Fig. 2.9, are illustrated for t = 1000. These trajectories are plotted on top of the nal locations of the other 245 cells at time t = 1000. Note that the distance travelled by the particles is longer in both models when the local interaction distance is smaller. We further explore properties of Model A by considering multiple runs with the same initial conditions and setup. In Fig. 2.12, we consider six di erent runs, all with an interaction distance D = 10 pixels and N = 250 cells at t = 1000 time steps. Note that in comparing the six images the resulting distribution of cells and the size and number of aggregations are similar, but there do not appear to be any stable aggregations which routinely form in the same location. In fact, the aggregations which form move over time as well as dissociate and form new aggregations. 32 D = 5 D = 10 (a) (b) M ode l A (c) (d) M ode l B Figure 2.9: (a), (b) Model A; (c), (d) Model B. Periodic boundary conditions; t = 1000; N = 250 particles. The interaction distance is D = 5 for (a), (c) and D = 10 for (b), (d). The initial conditions and parameters a, b are identical to those in Fig. 2.7. 33 # of neighbors within 5 pixels # of neighbors within 10 pixels (a) (b) M ode l A , D = 5 (c) (d) M ode l B , D = 5 Figure 2.10: Histogram illustrating how many cells have 0 to 10 neighbors within 5 or 10 pixels. Sub gures (a), (b) correspond to the simulation of Model A with D = 5 in Fig. 2.9(a). Sub gures (c), (d) correspond to the simulation of Model B with D = 5 in Fig. 2.9(c). In (a) and (c), we graph the number of neighbors that are within 5 pixels of each cell. In (b) and (d), we graph the number of neighbors that are within 10 pixels of each cell. 34 D = 5 D = 10 (a) (b) M ode l A (c) (d) M ode l B Figure 2.11: The same ve cells traced over 1000 time steps for (a), (b) Model A, and (c), (d) Model B. Periodic boundary conditions; t = 1000; N = 250 particles. The interaction distance is D = 5 for (a), (c) and D = 10 for (b), (d). The initial conditions and parameters a, b are identical to those in Fig. 2.7. Model A allows for more persistence in directed movement. 35 (a) (b) (c) (d) (e) (f) Figure 2.12: (a){(f) Six realizations of Model A with the same initial conditions and parameters a and b, using D = 10 and N = 250 cells at time t = 1000 time steps. While precise locations of aggregations vary, the aggregation sizes and numbers are similar. Varying the interaction distance D, we can gain further insight about the properties of this model. In Fig. 2.13, we consider three di erent interaction dis- tances D =10, 20, and 40 at t = 1000 time steps with the same initial conditions but di erent boundary conditions. In Fig. 2.13(a){(c), cells are constrained to a circle, whereas in Fig 2.13(d){(f), cells are simply restricted to R2 (note that in these cases not all cells are necessarily visible in the viewing window). In comparing the results for D = 10 and D = 20, we observe that the distribution and number and size of aggregations of cells is very similar. For D = 40, even the locations of the aggregations is almost identical. This is typical for simulations where D = 40; however, it is also possible for three aggregations of cells to form. The positions of the four aggregations and the distances between aggregations are generally slightly 36 di erent. In both cases, all cells for D = 40 typically end up in an aggregation. (a) (b) (c) (d) (e) (f) Figure 2.13: A study on the e ect of interaction distance D. We consider (a), (d) D = 10, (b), (e) D = 20, and (c), (f) D = 40 with the same initial conditions and same parameter a, b and N = 250 cells at time t=1000. For (a)-(c), the cells are constrained to a 200 pixel diameter circle. For (d)-(f), the cells are allowed to move anywhere in the R2 plane. As D increases, we observe fewer aggregations but with more particles. 2.5 Conclusions In this chapter, we sought to answer a series of four questions regarding the seemingly random motion of and communication between cyanobacteria. In our at- tempt to answer these questions, we presented two mathematical models for study- ing the local interactions during the motion of Synechocystis sp. One model (B) assumed that a cell chooses to move in random directions and another model (A) assumed that a cell may choose to move in the direction of ones of its neighbors. Comparing these models, we can discuss question 1. While we were not able to 37 completely isolate whether the chosen direction of cellular movement is completely random or not, in model A, we do observe aggregations and interactions slightly more consistent with experiments. Additionally, our assumptions of model A o er a plausible mechanism by which aggregations form, a possible answer to question 4. To approach questions 2 and 3, for both models we assumed that cells only move if a cell has a neighbor within some interaction distance D. We brie y compared this to a model of completely random motion as well as varied the interaction distance D to study the e ect of such a distance. The movement observed with models A and B was much more consistent with experimentally observed cell movement than with the purely random model. However, it is still possible that the interaction distance D should instead be modeled as an independent property of each cell, possibly even varying with time. Beyond these questions, there are similar characteristics of the model simula- tions and experimentally observed cell movement. True to the experimental obser- vations, both models never reach a stationary steady state. The cells, unless stuck to the boundary with no neighbors within interaction distance D, are in a constant state of motion. Cells that do not permanently adhere to the boundary are capable of moving again at a later time. Both models include the persistence in the motion; yet, it is the possible that the direction of motion will change at any time. While the intercellular signaling in Synechocystis sp. is not well understood, it is thought to be present. The local interaction model we have developed does not specify how the signaling takes place. It is possible that signaling depends on the proximity or connectedness of pili of individual cells. Additionally, cells might 38 release compounds which signal to other cells which direction to move. For instance, it has been shown that cAMP may be part of the signaling pathway required for motility [12]. The interaction distance in our model could be related to the time scale for di usion of such a molecule, or possibly be related to the length of a pilus. If the latter is true, our model suggests that cells that have been genetically modi ed to have longer pili might form larger aggregates of cells, a possible biological implication to answer question 3. Both models, the local-interactions model and the neighbor-dependent random model, produce results that better capture the experimental results than a model of a simple random motion. These models provide a possible explanation of the ob- served quasi-random behavior. The interaction distances for simulations considered in this paper seem to provide upper and lower bounds for reasonable behavior. Fur- ther study of this parameter, in connection with density and experimental results, may provide more insight into the relevant communication length-scale between cells. Understanding the local interactions between cells may turn out to be the key ingredient in understanding the mechanisms that control phototaxis. This might ul- timately lead to understanding the source of the observed dynamical patterns (such as the nger formation, delayed motion, etc.). We have carried out a preliminary study of the integration of local-interaction models with global forcing (due to the light source) in [32], which we discuss further in the next chapter. 39 Chapter 3 Phototaxis: incorporating global forcing due to the light source 3.1 Introduction In this chapter, we extend our model of local interactions from Chapter 2 to include global forcing due to light. This is done by studying an activation process of cells as they become a ected by the presence of light. Several published mathematical models of phototaxis focused on light, in ef- forts to replicate the experimentally observed emerging dynamical patterns of cells. Burriesci and Bhaya derived an agent-based model where cells follow a biased ran- dom walk in the direction of light with preferential movement toward regions previ- ously occupied by cells [15]. They also analyzed experimental data to determine the average velocity of cells in various phases of their motion, i.e. at the front of a nger, back of a nger, or along an edge, and determined that at least one cell is motile in 100% of aggregates of 8 or more cells, about 70% for groups of 4 to 6 cells, and less than 10% for groups of 3 or fewer cells. The model in [15] did not incorporate group size or any form of quorum sensing. All cells ultimately move in the direction of light, producing nger-like projections. Levy and Requeijo derived a cellular au- tomaton model and a related system of stochastic di erential equations that model the position of particles, the surface on which the particles move, and the excitation level of particles, a novel concept where cells become excited based on the average 40 excitation of their neighbors [55]. The excitation level determined the ability of a cell to sense light. Once the excitation of cells passes a critical threshold, their motion is directed towards the light source with the addition of a two-dimensional Brownian motion. Levy and Requeijo also were able to derive the continuum limit in the form of a system of partial di erential equations from their stochastic particle system [54]. Ha and Levy extended the Cucker-Smale ocking model to produce a particle system accounting for cell position, velocity and excitation [36]. With their particle model, they are able to derive a kinetic model, from which they are able to derive a Vlasov-McKean equation coupled to a reaction di usion equation. They proved that ocking behavior occurs for the particle model as well as for the kinetic model (under the assumption that all particles are fully excited). While the previous works capture the formation of nger-like projections, they do not capture the observed quasi-random motion of cells. In fact, these models typ- ically produce simulations where all cells ultimately move toward light or remain stationary, contrary to what is experimentally observed. Typically cells do not all ock together, and some cells become stationary singlets or align along edges. There may be various reasons why cells do not to migrate towards light; for instance, it is possible that not all cells are capable of identifying the direction of light. Alter- natively, it is possible that cells might \prefer" to move in di erent directions than the direction of light, due to signaling from neighboring cells, e.g. In this chapter, we extend our mathematical model of local interactions from Chapter 2 to include global forcing due to a light source. Instead of focusing on nger formation, we examine the emerging behavior in the entire population. 41 This chapter is organized as follows. We begin in Section 3.2 by extending the model of Chapter 2 to include the possibility of a cell sensing a global force due to light. We then hypothesize that not all cells are able to identify the direction of the light source and rewrite our model accordingly. Simulations of these global forcing models are presented in Section 3.3; these results are published in a conference proceedings [32]. In Section 3.4, we develop a new model of the activation process by which cells become excited and develop the ability to preferentially move in the direction of light. The results obtained with this model are provided in Section 3.5. Finally, in Section 3.6 we discuss the biological and mathematical implications of our models. 3.2 Global forcing model for phototaxis We start by overviewing several experimental observations that supplement what we presented in Section 2.2. 1. Finger formation and global forcing due to light. The dominant char- acteristic of phototaxis is the formation of nger-like projections. These pro- jections typically, for phototactic cells, reach out in the direction of the light source. In this way, light can be thought of acting as a global force on the system. Examples of these nger-like projections can be seen in Figure 2.1. 2. Activation process. Before any cell begins moving in the direction of light, it undergoes an activation process. Cells are observed to form aggregations. When these aggregations reach critical sizes, cells are able to sense the direction 42 of light and move in that direction. For small aggregations of cells, little movement is observed [15]. 3. Inhomogeneity. We observe that not all cells respond identically to the light source. While we assume that the cells under consideration, Synechocystis sp., all have the same genome (with minor mutations), it is possible that they display di erent phenotypes, having di erent ages or light-sensing abilities. This might be a possible explanation to why not all cells necessarily move in the direction of light, as we demonstrated in Chapter 2. 3.2.1 Model assumptions We begin by making the same assumptions as we did for our model of local interactions (the so-called Model A in Chapter 2). We assume that cells can identify the locations of their neighbors, but we do not assume anything about underlying communication mechanisms between cells { possibly by signaling, local force sensing, or physical interaction. For the purpose of this model, the actual mechanism of communication is of no interest. Again, the size of the neighborhood, in which cells can communicate with each other, will be taken as one of the model parameters. Additionally, we assume that cells are not necessarily all alike, i.e., it is possible for cells to have di erent traits. It is experimentally observed that not all cells move in the direction of the light, even after long periods of time and especially when cells are in isolation. This implies the possibility that not all cells are capable of identifying the direction of the light. 43 Global forcing can be incorporated into our stochastic model of local interac- tions in di erent ways. Accordingly, we propose two models that incorporate the individual response of cells to the light. In the rst model (Model H) we assume that all cells are capable of identifying the direction of the light; that is, the cells are homogeneous in their directional light sensitivity. In the second model (Model I) we assume that only a subset of the colony is capable of identifying the direction of light; that is, the cells are inhomogeneous in their directional light sensitivity. In both models, cells that are capable of sensing the direction of light are not assumed to automatically move in that direction. Rather, when moving, in addition to the local-interactions rules, they have the additional option of moving towards light with a xed small probability. We illustrate the four possible cellular motions in Figure 3.1. (a) (b) (c) (d) Figure 3.1: (a){(c) Local interaction model and (d) Global forcing in direction of light. Model particles can exhibit (a) persistence, (b) stationary behavior, (c) movement toward neighboring cells, and (d) movement in the direction of a light source. 44 3.2.2 Model H: Homogeneous response to directional light force In this model, every cell is capable of sensing the direction of light. We modify equation (2.1) to allow the cells to move in the direction of light, 0, with probabil- ity c. Additionally, each cell can choose to persist in its previous direction of motion with probability a. Cells can also change their velocity from stationary to moving, or moving to stationary, with probability 1 b. To formulate the model, we consider N particles in R2. Each particle i has a direction i(t) 2 ( ; ] and a velocity vi(t) 2 f0; 1g. The time is assumed to be discrete with t 2 N. The position of particle i at time t is denoted by xi(t). Similar variables are considered in a recent discrete model of chemotaxis [63]. We update the direction of each particle i according to i(t+ 1) = 8 >>>>>>>>< >>>>>>>>: i(t); with probability a, 0; with probability c, j; with probability 1 a cni for j s.t. xj 2 BD(xi), (3.1) where BD(xi) is a ball of radius D centered at xi, ni is the number of particles with position xj contained in BD(xi), i.e. jjxi(t) xk(t)jj2 < D, and j is the angle of the vector pointing from particle i to particle j. When a particle chooses to move in the direction of a neighboring particle, the probability of choosing any neighboring particle is equal. The velocities are updated according to the following rule: vi does not change with probability b. With probability 1 b the velocity vi switches 45 between 0 and 1. This is formulated as vi(t+ 1) = 8 >>>< >>>: vi(t); with probability b, 1 vi(t); with probability 1 b. (3.2) In this way, the velocity is either 1 (in which case the cell continues to move) or 0 (in which case it stays in place). After calculating the angle and the velocity for every particle, the positions fxig of the cells are updated according to xi(t+ 1) = xi(t) + vi(t) 2 6 6 4 cos( i(t+ 1)) sin( i(t+ 1)) 3 7 7 5 : (3.3) Note that in this case, we assume that the light source is far away so that the direction of light is identical for all cells on the plate. 3.2.3 Model I: Inhomogeneous response to directional light force In the model of inhomogeneous light sensitivity, only a subset of particles are capable of sensing the direction of the light source. To model this, particles are randomly assigned to be of type leader or type normal. A particle of the leader subtype updates its angle of motion according to equation (3.1). A particle of the normal subtype updates its angle of motion according to the rules of motion from the local interaction Model A in equation (2.1). Note that by setting c = 0, equation (3.1) reduces to equation (2.1). The ratio of leader to normal cells is a model parameter. Its e ect is studied in the numerical simulations in the next 46 section. 3.3 Simulation results of the models of Section 3.2 In this section we present results of simulations of our global forcing models. In all cases, we assume that the cells are constrained to a 200 pixel diameter disc in R2. Cells are allowed to move in any direction that is allowed by the model within the disc but are not allowed to cross its boundary to the exterior of the domain. Such a constraint is consistent with the observed behavior of cells on the wetted area of a plate on small time-scales. On larger time-scales, in the absence of a directed light source, cells di use radially. In each simulation, we set the number of particles N to be either 250 or 500. The particles are represented in each simulation by circles of radius 2 pixels. For each set of parameters, we consider two di erent interaction distances: D = 5 and D = 10 pixels. The interaction distance is shown with a red bar of the appropriate pixel length in the upper left corner of each gure. In all simulations we assume that the probability of directional persistence, a, is 0.8. The probability of switching velocity from 1 to 0 or 0 to 1, 1 b is set as 0.5. For both global interaction models H and I, the probability of a cell moving in the direction of light at any time step, c, is 0.005. For the global interaction Model I, unless otherwise stated, the probability of a cell being able to sense the direction of light is 0.5, i.e., the ratio of leader to normal cells is 1:1. The group of leader cells are identi ed as un lled circles in the simulations. For convenient reference and comparison, we again display one gure with 47 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 3.2: Local interaction Model A at t = 0; 50; 350 time steps with (a)-(c) D = 5, N = 250, (d)-(f) D = 10, N = 250, (g)-(i) D = 5, N = 500, and (j)-(l) D = 10, N = 500. 48 results from the local interaction Model A from Chapter 2. In the local interaction model, cells are capable of moving toward a random neighbor, persisting in their chosen direction, or being stationary. Simulations of the local interactions base model for N = 250 and 500 and for D = 5 and 10 are shown in Figure 3.2. The snapshots for each simulation are taken at times t = 0; 50; and 350. The initial positions, directions, and velocities of the cells are identical for Fig. 3.2(a) and Fig. 3.2(d) with N = 250 and for Fig. 3.2(g) and Fig. 3.2(j) with N = 500. The same initial conditions are used in all of the following simulations. Note that in the case D = 5, cells become temporarily xed on the edge of the circle, visible in Fig. 3.2(b), Fig. 3.2(c), Fig. 3.2(h) and Fig. 3.2(i). These cells occasionally move back into the domain, but generally, spend most of the time on the boundary. We also observe large aggregations forming for D = 10 in Fig. 3.2(f) and Fig. 3.2(l) and small aggregations forming for D = 5 in Fig. 3.2(c) and Fig. 3.2(i). All images illustrate the observed quasi-random motion of cells that is observed experimentally. In the global Model H, all cells are aware of the direction of light but only have a 0.5% chance of choosing to orient in the direction of the light at each time step. Note that if a particle is moving in the direction of light (or in any other direction), with each progressive time step, the cell will persist in that direction with probability 0.8. Simulations of this joint global-local interaction model are shown in Figure 3.3. The initial conditions fxi(0); i(0); vi(0)g are the same for each set of values of N and D as in Figure 3.2 (see Figs. 3.3(a), 3.3(d), 3.3(g) and 3.3(j)). The direction of light in these simulations is at the lower right corner ( = =4). For D = 5, we still see cells that become xed on the edge of the circle, even for 49 t = 1000 (see Figs. 3.3(b), 3.3(c), 3.3(h) and 3.3(i)). For D = 10, we still observe the formation of large aggregations. In comparing Fig. 3.3(c) with Fig. 3.3(f) and Fig. 3.3(i) with Fig. 3.3(l) at t = 1000, the mass movement in the direction of light is clearly observable for D = 5 but not in the case of a larger neighborhood, i.e., D = 10. The migration toward light appears to be greater for the smaller number of particles N = 250 compared with N = 500 as can be seen in Figs. 3.3(c) and 3.3(i). The global Model I combines the local interaction model with only a fraction of cells that are allowed to intentionally pursue the direction of the light source as in global Model H. Simulations of this global model with leader to normal cell ratio of 1:1 are shown in Figure 3.4. The initial conditions of variables fxi(0); i(0); vi(0)g are identical to those in Figs. 3.4(a), 3.4(d), 3.4(g), and 3.4(j). The cells which are aware of the direction of light are identi ed in the simulation as open red circles. In these simulations, as the number of time steps increases, one expects to see general migration of cells toward the lower right hand quadrant of each circle as happened with H. Similar to Figure 3.3, this trend is not apparent at t = 500, in Figs. 3.4(b), 3.4(e), 3.4(h) and 3.4(k), or at t = 1000 and D = 10, in Figures 3.4(f) and 3.4(l). It can be observed in the distribution of red leader cells for D = 5 at t = 1000 in Fig. 3.4(c) and Fig. 3.4(i). The distribution of normal cells appears to be uniform in each of these snapshots. Large time simulations of each model are shown in Figure 3.5. Note that for the local interaction model alone, in Figs. 3.5(a), 3.5(d), 3.5(g) and 3.5(j), even though cells are still moving, the pattern of cells does not change signi cantly between t = 3000 and t = 9999. In the global Model H (shown in Figs. 3.5(b) and 3.5(e) 50 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 3.3: Global forcing Model H at t = 0; 500; 1000 time steps with (a)-(c) D = 5, N = 250, (d)-(f) D = 10, N = 250, (g)-(i) D = 5, N = 500, and (j)-(l) D = 10, N = 500. The direction of light is = =4. 51 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 3.4: Global forcing Model I with leader (red ) to normal (black ) cell ratio 1:1 at t = 0; 500; 1000 time steps with conditions (a)-(c) D = 5, N = 250, (d)-(f) D = 10, N = 250, (g)-(i) D = 5, N = 500, and (j)-(l) D = 10, N = 500. The direction of light is = =4. 52 with D = 5), while most of the cells aggregate along the side arc of the circle that is closest to the light source, cells are still capable of escaping back into the domain. Whereas, in Fig. 3.5(h) and Fig. 3.5(k) with D = 10, most of the cells aggregate along the leading edge as the number of time steps increases. In the global Model I with D = 5 (shown in Figs. 3.5(c) and 3.5(f)), the leader-type cells move towards the leading edge but do not appear to gather along the edge in large numbers. Instead, they end up spreading throughout the domain. The normal-type cells are spread somewhat uniformly. These patterns are indicative of the observed quasi-random motion of cells. For D = 10 at t = 3000, the aggregations of leader- and normal- type cells are somewhat uniform in Fig. 3.5(i), but most aggregates tend towards the light source as t increases as can be seen in Fig. 3.5(l). Finally, we study the impact of the ratio of leader to normal cells in the global Model I. Simulations are shown in Figure 3.6. We consider three cases in which the ratios are 3:1, 1:1, and 1:3. In addition, we consider the case where all cells are leader cells, as in global Model H. The initial conditions fxi(0); i(0); vi(0)g are taken to be the same as before. They are shown in Figs. 3.6(a), 3.6(d), 3.6(g), and 3.6(j). In all simulations, D = 5 and N = 250. The results of the simulations are consistent with the characterization of the quasi-random motion. We also observe that cells get stuck along the entire edge of the circle. A migration of leader cells toward light, in all three ratios, (shown in Figs. 3.6(b), 3.6(e), and 3.6(h)) is somewhat visible and is very apparent in Figs. 3.6(c), 3.6(f), and 3.6(i). In each of these cases, the distribution of normal cells appears to be uniform. Comparisons of these simulations with experimental data to the simulation where all cells are leader-type 53 Local model A Global model H Global model I (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 3.5: Long time simulations (t = 3000 and 9999) of each model for N = 250 particles with (a)-(f) D = 5 and (g)-(l) D = 10. The direction of light is = =4. 54 cells supports the hypothesis that cells may not be all alike. It seems feasible that only a portion of the cells are capable of identifying the direction of the light source. 3.4 An activation model for phototaxis In Section 3.3, we illustrated the possible implications of assuming an inho- mogeneous light sensitivity. In this section, we consider incorporating an activation process in which particles can self-propel their motion or follow other particles based on an internal excitation process. During the activation process, cells become ex- cited by the presence of light and aggregate prior to moving in the direction of light. Aggregations seem to enhance the apparent excitement level and sensitivity of a cell to light, as indicated by the data analysis in [15]. Following [55], we model this process by assigning each particle i an internal excitation, denoted Si(t), at each time t. The criteria we follow when constructing Si(t) are di erent than in [55]. We construct the excitation process as a positive value that meets the following criteria: 1. All particles are initiated in an unexcited normal state, where the internal exci- tation is below a critical threshold, K. Above this critical threshold, particles become \excited" and switch to the leader subtype. 2. Excitation of particles builds in the presence of light and decays in the absence of light. 3. The excitation values must saturate. Furthermore, the excitation level of single isolated particles, with no detectable neighbors, should saturate below the 55 (a) 1 leader : 3 normal (b) 1 leader : 3 normal (c) 1 leader : 3 normal (d) 1 leader : 1 normal (e) 1 leader : 1 normal (f) 1 leader : 1 normal (g) 3 leader : 1 normal (h) 3 leader : 1 normal (i) 3 leader : 1 normal (j) All leader type (k) All leader type (l) All leader type Figure 3.6: Simulations of global Model I at t = 0; 1000; and 3000 time steps for D = 5 with various ratios of leader (red ) to normal (black ) cells: (a)-(c) 1 leader : 3 normal, (d)-(f) 1 leader : 1 normal, (g)-(i) 3 leader : 1 normal and (j)-(l) All leader cells, i.e. global Model H. The direction of light is = =4. 56 critical threshold. The excitation of particles in aggregations, with detectable neighbors, should saturate at a level above the critical threshold. Previous models of activation were developed such that the excitation of a given particle is driven towards the mean excitation of neighboring particles. In- stead of being able to determine the mean excitation of surrounding particles, we assume that a particle can determine the total excitation of surrounding particles and compare this sum to its internal excitation. To account for this, we de ne a quantity i as i = P j;xj2BD(xi) Sj Si ; (3.4) where BD(xi) is a disc of radius D of centered at the location of particle i, xi. Note that this disc always contains xi and, hence, i is always greater than or equal to one. We can interpret i to be the e ective number of neighbors that a particle i senses relative to the excitement levels of surrounding particles and the internal excitement of particle i. We describe the activation process by an ordinary di erential equation. Sup- pose the maximum constant rate of excitation production is L and that this is attenuated based on the e ective number of neighbors of a particle i and based on the total excitation level within detection distance D. The total excitation surround- ing a particle is calculated by iSi. To attenuate this maximum rate, we multiply by Michaelis-Menten terms in i and iSi, with Michaelis constants K1 and K2, respectively. By de nition, the Michaelis constant is the value at which the system attains half of the maximum rate. We also allow the excitation Si to decay at a rate 57 proportional to Si S0, with rate constant , where S0 is a baseline, below which, the excitation increases and above which the excitation decreases. The resulting ODE for the excitation is d dt Si = L i K1 + i iSi K2 + iSi (Si S0): (3.5) Since accuracy is of no concern, we use a simple forward Euler solver to ap- proximate the solution of (3.5), allowing the unit time step to match our increments used in other simulations. Accordingly, the expression we evaluate is Si(t+ 1) = Si(t) + L i(t) K1 + i(t) i(t)Si(t) K2 + i(t)Si(t) (Si(t) S0): (3.6) In the simulations shown in Section 3.5, we x a few of the parameter values and explore reasonable ranges for the other parameters. We set the production rate at L = 0:1 and the decay rate as = 0:02. The ratio of these parameters, as well as the baseline S0, determines the maximum attainable value of Si. We x the baseline at S0 = 1. In this way, the maximum saturation value of Si for all particles is 6. In the absence of light, we can set L = 0, but in the simulations of Section 3.5, we always assume that a light source is present. We initialize the excitation values at Si(0) = 1:1. For the Michaelis constant K1, we consider the values 2, 4, and 8. The value of K1 represents the ability of particles to move in the direction of light based on the e ective number of neighbors. We justify our chosen range of values by noting that, experimentally, cells are observed to have a high probability 58 of moving in the direction of light after aggregating in groups of size 6 or larger and a low probability of moving in the direction of light if it is in a small aggregation of 3 or fewer cells [15]. We x the second Michaelis constant K2 = 10; in this way, if a particle is surrounded by multiple excited particles, it may attain the maximum rate of excitation production. Being surrounded by a few unexcited particles will reduce the rate of excitation production. For the excitation threshold we consider the values K = 1:5 and 2.5. Note that we do not assign units to excitation. While there are some biological systems for which excitation is well-established, in our context its existence is speculative, and therefore it does not correspond directly to any physical or measurable quantity. 3.5 Simulation results of activation model We perform multiple simulations using our activation model (3.5). The initial conditions used in all simulations of our activation model are identical to those used in Figure 2.5 and throughout Chapter 2, with the additional constraint that the initial excitation Si(0) = 1:1 for all i. Furthermore, we x the persistence probability a = 0:8, the stopping probability b = 0:5, and the light detection probability, when applicable, c = 0:01. We compare the e ect of varying the interaction distance and consider the cases of D = 5 and D = 10. The length scale of the interaction distance D is shown in the upper left corner of each picture around a stationary cell (which is part of the legend and not part of the simulation). In each simulation, as before, the unexcited normal particles are shown as black dots and the excited leader particles 59 are illustrated with red open circles. The direction of light in each simulation is towards the bottom left-hand corner of each image; that is, 0 = 3 =4. In the rst simulation of our activation model, we consider the evolution of our particle system for the various values of the e ective neighbor number Michealis constant K1 = 2, 4, and 8. In Figure 3.7, we simulate 250 particles, constrained to a circle of diameter 200, with a neighbor detection distance D = 5 and an excitation cuto K = 1:5 at times t = 50, 350, and 1000. We observe that for K1 = 2 in Fig. 3.7(a){(c), approximately half of the particles become excited (the particles shown in red circles) and move in the direction of light. For K1 = 4 in Fig. 3.7(d){ (f), less than one quarter of particles become excited and move in the direction of light. For K1 = 8 in Fig. 3.7(g){(i), very few particles become excited; this is most likely due to the relatively low number of particles present in the small aggregations. Next, we consider the long-time behavior of our particle system for the various values of the e ective neighbor number Michealis constant K1 = 2, 4, and 8, two excitation cuto s K = 1:5 and 2:5, and two interaction distances D = 5 and 10. In Figure 3.8, we simulate 250 particles, constrained to a circle of diameter 200 at time step t = 3000. For the lower cuto K = 1:5 in Figures 3.8(a){(f), we observe that particles are able to aggregate and move in the direction of light in all but one simulation. For the larger neighbor detection distance D = 10, more particles form aggregations, resulting in more particles becoming leaders. For the higher cuto K = 2:5 in Figures 3.8(g)-(l), particles are able to aggregate and activate for D = 10 in (j){(l), but the threshold K = 2:5 is clearly too high for particles to detect the light in any reasonable manner. 60 (a) (b) (c) (d) (e) (f) (g) (h) (i) Figure 3.7: The activation model with (a){(c) K1 = 2, (d){(f) K1 = 4, and (g){(i) K1 = 8 for excitation cuto K = 1:5. N = 250 particles are constrained to a circle of diameter 200. Shown are simulation results at times t = 50; 350; 1000. Identical initial conditions are used in all simulations as well as parameters D = 5, a = 0:8, c = 0:01 and b = 0:5. The direction of light is = 3 =4. 61 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 3.8: The activation model with at long-time t = 3000 for (a){(f) cuto K = 1:5 and (g){(l) K = 2:5. The detection distance is D = 5 in (a){(c) and (g){(i) and D = 10 in (d){(f) and (j){(l). The rst column has K1 = 2, the second column has K1 = 4, and third column has K1 = 8. In each simulation, N = 250 particles are constrained to a circle of diameter 200. Identical initial conditions are used in all simulations as well as parameters a = 0:8, c = 0:01 and b = 0:5. The direction of light is = 3 =4. 62 In Figure 3.9, we illustrate that the excitation values Si(t) are in fact bounded. Each gure shows the time dynamics of the excitation of all particles. The parame- ters and locations of each image correspond to the plots in Fig. 3.8. As you can see, for systems where many particles become excited, here D = 10 in Fig. 3.9(d){(f) and (j){(l), the excitation levels do in fact saturate at a value below 6.5, the expected maximum. Note that they saturate below this level because the Michaelis-Menten terms reduce the excitation production rate. Furthermore, for parameter sets where relatively few particles become excited, the particles in the system do not reach the saturation level; instead, we observe oscillations of particles alternating between excited and unexcited states. Finally, we study the trajectories of individual particles. In Figure 3.10, we illustrate the trajectories of 5 particles over 1000 time steps for various parameter values. We only consider the excitation threshold K = 1:5. In a few of these images we are able to observe trajectories of particles moving in the direction of light as part of an aggregation of particles; for instance, see the green and pink trajectories of (a) and the pink and blue trajectories of (d) and (e). These patterns of motion are coherent with the experimental observations. 3.6 Conclusions and future directions In this chapter we extended the discrete-time, stochastically interacting par- ticle model from Chapter 2, which described local interactions between phototactic cells in which cells can stochastically choose a neighboring cell to follow, to include 63 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 3.9: Scatter plot of the excitation levels at long-time t = 3000 for (a){(f) cuto K = 1:5 and (g){(l) K = 2:5. The y-axis is excitation value Si(t) and the x-axis is time t. In each gure, the excitations for all 250 cells are plotted. The detection distance is D = 5 in (a){(c) and (g){(i) and D = 10 in (d){(f) and (j){(l). The rst column has K1 = 2, the second column has K1 = 4, and third column has K1 = 8. In each simulation, N = 250 particles are constrained to a circle of diameter 200. Identical initial conditions are used in all simulations as well as parameters a = 0:8, c = 0:01 and b = 0:5. The direction of light is = 3 =4. 64 (a) (b) (c) (d) (e) (f) Figure 3.10: Examples of cell trajectories over t = 1000 time steps for the activation model. The detection distance is D = 5 in (a){(c) and D = 10 in (d){(f). The rst column has K1 = 2, the second column has K1 = 4, and third column has K1 = 8. In each simulation, N = 250 particles are constrained to a circle of diameter 200. Identical initial conditions are used in all simulations as well as parameters K = 1:5, a = 0:8, c = 0:01 and b = 0:5. The direction of light is = 3 =4. 65 global forcing due to light. We varied several factors in each simulation: the inter- action distance, D, the number of particles in simulation, N , and the fraction of cells which were sensitive to moving in the direction of light. We also developed a mathematical model of the activation process, by which cells aggregate and become capable of moving in the direction of light. In simulations of this model we vary the model parameters for the interaction distance D, a Michaelis constant for e ective aggregate size K1, and the excitation cuto K, a threshold for determining whether or not a particle is capable of sensing the direction of light. We incorporated global forcing due to a light source into our local interaction model from Chapter 2 by stochastically allowing cells to move in the direction of light (in addition to the other local rules of motion). We considered two scenarios: one in which all cells are capable of identifying the direction of light and a second case in which only a fraction of the cells are capable of identifying the direction of light. When all cells are capable of identifying the direction of light, they all tend to migrate in that direction. However, on large time scales, for a smaller interaction distance it is still possible to see some quasi-random motion. This is consistent with the experimental data. When comparing all models on large time scales, a quasi-random motion seems to be more present in the models where not all particles are capable of detecting the direction of light. This issue should be further investigated experimentally, but the evidence presented in this chapter strongly suggests that there may be a substantial subpopulation of cells that is not capable of detecting the direction of the light source. 66 Our activation model, under certain parameter sets, is able to replicate the experimentally observed process of activation, where cells aggregate before moving in the direction of light. Furthermore, the internal excitation of each model particle saturates, a property which was not attained by previously developed models. With a relatively simple model using Michaelis-Menten kinetics, we are able to vary the relative size of aggregates in which particles must reside in order to identify and move in the direction of light. Further controlled experimental and numerical studies should be performed in order to determine what parameter set best describes the motion of phototaxis in Synechocystis sp. It is likely that such a parameter set would vary in time as cells display di erent motion characteristics in di erent regions. 67 Chapter 4 One-dimensional local interaction models 4.1 Introduction In this chapter, we consider a one-dimensional model of the local interactions model presented in Chapter 2. One-dimensional models are often used to simplify systems and make problems more tractable. Simulations of one-dimensional models are typically less computationally intensive. In addition, analytical techniques are more accessible in one-dimensional problems. Throughout the chapter, we follow the ideas of Baker, Yates and Erban who developed a microscopic model for studying cell migration on growing domains in one dimension [5]. Starting from a set of basic rules of motion, they develop a reaction di usion master equation (RDME), from which they derive a system of ordinary di erential equations (ODEs) describing the cell populations along an in- dexed number line. Finally a partial di erential equation is obtained as a limit of the ODEs. Baker et al. focused on a variety of signaling patterns and boundary conditions to study cell migration on growing domains. Following the methodology in Baker et al. [5], we derive a system of ODEs from a one-dimensional version of our stochastic discrete model from Chapter 2. We compare simulations of our stochastic 1-D particle system and the new system of ODEs to illustrate the strengths and weaknesses of each model. We also utilize 68 the system of ODEs to explore the parameter space, which we can now do more thor- oughly because solving the system of ODEs requires less computational resources than solving the corresponding stochastic particle system. 4.2 Local interaction model for two populations in one dimension We consider a set of particles which are constrained to a one-dimesional array. In order for this model to resemble our previous two-dimensional models, we allow particles to sense neighbors within a given distance and to choose a direction for their motion that is based on their neighbors. In the one-dimensional case, our model becomes a system of right-moving and left-moving particles on a line, as particles can only move to the right or to the left of their current positions. At every time step, each particle can move one bin to the left or to the right. We will later allow the particles to remain stationary. We begin our model development by looking at a set of particles, each with a given initial position and initial preferred direction. First, let us consider the e ective transitional probabilities for particles moving between bins indexed along a number line. Suppose that right-moving particles, that is, those that came to be in their current position by moving to the right, persist in their motion and continue to move to the right with probability a. With probability 1 a, a particle may choose a new direction of motion (which can be either to the right or to the left). The likelihood of a right-moving particle i at position xi moving left is weighted by the number of detectable neighboring particles located to the left, normalized by the total number of particles in both directions. Assume that 69 particles can detect particles located within distance D to the left or to the right. We let i denote the number of particles located within the neighbor detection distance D to the left of particle i and +i denote the number of particles within the detection distance D to the right of particle i. We illustrate this setup in Figure 4.1. Accordingly, the probability of a right-moving particle choosing to move left is (1 a) i =( + i + i ) and the probability of a right-moving particle choosing to move right is (1 a) +i =( + i + i ). Similarly, the likelihood of a left-moving particle i at position xi moving right is weighted by the number of detectable neighbor cells located to the right +i , divided by the total number of its neighboring particles moving in either direction ( i + + i ). The probability of a left-moving particle persisting to the left at the next time step is a. The probability of a left-moving particle choosing to move left is (1 a) i =( + i + i ). We note that we di erentiate between a particle persisting to move in a direction and a particle choosing in which direction to move (which can still be the old direction of motion). By adjusting the probabilities, we could limit the particles to persist in their direction or choose a new direction. Also note that the probability of choosing to move left and the probability of choosing to move right depend only on the location of neighboring particles and not on the previous direction of motion of the particle. xi !i"! "## $### !i+! "## $### Figure 4.1: A number line illustrating the cell position xi of particle i with detectable neighbors being the total number of particles which fall in the regions indicated by +i and i . In this particular image, the neighbor detection distance D = 6. 70 To formulate a model, consider a set of N particles. The position of each particle i at time t is denoted by xi(t). We start by assuming that all particles are moving at a xed speed v, so that at every time step t the displacement is x = v t. In the next section, we will relax this condition. At every time- step, a particle can either persist in its previous direction of motion or choose a new direction; in the latter case the new direction can be either right or left, one of which must be its previous direction. When a particle moves to the right, its position is increased by x. When a particle steps to the left, its position is decreased by x. Suppose that the probabilities of choosing either direction are proportional to the number of neighbors in the D > 0 neighboring bins in either direction. Scaling t = 1, we have xi(t+ 1) = 8 >>>>>>>>>< >>>>>>>>>: xi(t) + [xi(t) xi(t 1)]; with probability a; xi(t) + x; with probability (1 a) +i +i + i ; xi(t) x; with probability (1 a) i +i + i : (4.1) where the number of particles to the right of bin xi that can be sensed by particle i is +i = DX m=1 NX j=1 xj ;xi+m x; (4.2) and the number of particles to the left of bin xi that can be sensed by particle i is i = DX m=1 NX j=1 xj ;xi m x: (4.3) 71 Here xi;xj is the Kronecker delta function; that is, xi;xj is 1 if xi = xj, and 0 otherwise. Note that this discrete stochastic system (4.1) is not Markovian; the position at time t+ 1 depends on the position at time t 1, not just the position at time t. For this setup, a variety of boundary conditions can be considered. In this chapter, we consider a xed number of particles on a xed interval with periodic boundary conditions. 4.2.1 Derivation of the Reaction-Di usion Master Equation (RDME) We begin by considering the number of particles in each bin, as opposed to the position of each particle; the former is more practical for systems with large numbers of particles. The number of particles ni(t) in bin i at time t is given by ni(t) = NX j=1 xj(t);i x: (4.4) Suppose that the particles are distributed between k bins. The number of particles in bins i = 1; : : : ; k is given by the vector: ~n = [n1; n2; : : : ; ni; : : : ; nk]: (4.5) Note that the vector ~n and each of its components ni are functions of time t, which has been suppressed here for the sake of brevity. At every time step, each particle has a preferred direction, in which it can 72 persist with probability a. Let us divide the particles into two populations: right- moving and left-moving. Let ri(t) denote the number of particles in bin i at time t that moved to the right to enter bin i. That is, ri(t) is the number of particles which were in bin i 1 at the previous time step t 1 and are in bin i at the current time t. Let li(t) denote the number of particles in bin i that moved to the left to enter bin i. These particles were in bin i + 1 at the previous time t 1. Observe that at each time step, the total number of particles in bin i is given by the sum of the right-moving particles and the left-moving particles: ni = ri + li: (4.6) We can now use eqs. (4.2), (4.3) and (4.4) to rewrite our expressions for the number of detectable neighboring particles in each direction: i = PD m=1 ni m and +i = PD m=1 ni+m. To simplify the notations, we de ne neighbor-weighted probabilities cri = (1 a) + i =( i + + i ) and c l i = (1 a) i =( i + + i ), which are the probabilities of a particle choosing to move to the right or the left, respectively. These notations are summarized in Figure 4.2. ri l i l i!1 ri+1 ci l a + ci l a + ci r ci r Figure 4.2: A sketch of the notation used for describing the motions of particles. Not shown in this gure are the arrows showing particles in bin i 1 moving into population ri and showing particles in bin i+ 1 moving into population li. a is the persistence probability and cri and c l i are the neighbor-weighted probabilities of a particle in bin i choosing to move to the right and left, respectively. 73 We de ne a joint probability density function, P (~r;~l; t), which describes the probability of the system attaining each state f~r;~l; tg for a given set of initial con- ditions. To describe the evolution of P (~r;~l; t) we must consider all ways in which particles can enter and leave the state f~r;~l; tg. This is done by summing over all of the bins i, accounting for methods by which a particle can leave a bin i. Motion of a particle in either direction results in a change in the state of the system. Depending on the present state of this system, this may a ect evolution into or out of the state f~r;~l; tg. We account for these items separately in Eq. (4.7) with items ? particles moving to the right to enter state f~r;~l; tg, ? particles moving to the right to leave state f~r;~l; tg, fi particles moving to the left to enter state f~r;~l; tg, and fl particles moving to the left to leave state f~r;~l; tg. Item is in place to account for the boundary conditions. @P @t (~r;~l; t) = k 1X i=1 h Probability a particle moves out of bin i to right to enter state f~r;~l; tg | {z } ? Probability a particle moves out of bin i to right to leave state f~r;~l; tg | {z } ? i + kX i=2 h Probability particle moves out of bin i to left to enter state f~r;~l; tg | {z } fi Probability a particle moves out of bin i to left to leave state f~r;~l; tg | {z } fl i + Boundary terms | {z } : (4.7) 74 To de ne each term in the RDME (4.7), we rst formulate particle-transposing operators which allow us to examine states that are adjacent to the state f~r;~l; tg. That is, only one particle is displaced from the state f~r;~l; tg. For this to happen, an additional particle must be present in some position i and absent from bin i 1 or i + 1. We denote these particle-transposing operators with J+i;r and J i;l; the operators are maps from Rk to Rk, operating on the distribution of particles ~r or ~l and producing a new distribution where the number of particles in bin i has increased by one and the number of particles in bin i + 1 or i 1, respectively, has decreased by one. As just stated, the particle-transposing operator for a right- moving particle moving to the right is J+i;r. The rst subscript indicates the position. The second subscript indicates the particle type that is moving. The superscript indicates the direction in which the particle is moving, right with + or left with . In equation (4.8) for operator J+i;r acting on ~r = [r1; : : : ; ri; : : : ; rk], there is an extra particle in bin i and one less particle in bin i+ 1 compared with ~r. J+i;r ~r = r1; : : : ; ri 1; ri + 1; ri+1 1; ri+2; : : : ; rk : (4.8) Similarly for the operator J i;l, which allows a left-moving particle in bin i to move to the left, there is an extra particle in bin i and one less particle in bin i 1 compared with ~l. This is expressed as J i;l ~l = l1; : : : ; li 2; li 1 1; li + 1; li+1; : : : ; lk : (4.9) 75 Note that the particle-transposing operators in eqs. (4.8) and (4.9) map from Rk to Rk and act only on vectors ~r and ~l. We need two more operators which allow particles to move between right-moving and left-moving populations. Such operators account for particles that change the direction of their motion. In equation (4.10), describing operator J+i;l, one additional left-moving particle is located in bin i so that it can move to the right into bin i+ 1. J+i;l ~r;~l = r1; : : : ; ri+1 1; : : : ; rk; l1; : : : ; li + 1; : : : ; lk : (4.10) To describe the operator J i;r, one additional right-moving particle is located in bin i and moves left to bin i 1. J i;r ~r;~l = r1; : : : ; ri + 1; : : : ; rk; l1; : : : ; li 1 1; : : : ; lk : (4.11) Note that particle-transposing operators J i;r and J + i;l act on [~r;~l], mapping from R 2k to R2k. With the above operators and probabilistic rates of transitioning to the right or left de ned, we are able to determine each component of the RDME (4.7). For ?, we want to account for the probability of a particle in bin i moving to the right to enter state f~r;~l; tg. The rate at which a right-moving particle moves right is a+ cri . Note that cri depends on the position of other particles and so we need to be a little more careful in handling this term. For a particle to move out of bin i to the right so that the new state is f~r;~l; tg, the particles must be arranged as in Figure 4.3. 76 Right-moving particles : : : ri 1 ri + 1 ri+1 1 ri+2 : : : Left-moving particles : : : li 1 li li+1 li+2 : : : Total : : : ni 1 ni + 1 ni+1 1 ni+2 : : : Figure 4.3: Setup of particles for J+i;r[~r] To calculate the neighbor-dependent weight which would allow a particle in bin i to move to the right, we sum over the neighbors to the right within neighbor detection distance D and normalize by dividing by the total number of neighbors within neighbor detection distance D. If we use the previous formula +i = PD j=1 ni+j with values from the current state, the number of neighbors to the right of bin i is exactly +i 1. The number of neighbors to the left, i , has not changed. Hence, for this setup, the weight for moving to the right is ( +i 1)=( + i + i 1). We can now calculate the probability of a right-moving particle moving to the right out of bin i. The probability that the particles are in the con guration shown in Fig. 4.3 is P (J+i;r[~r];~l; t). The number of particles which could move to the right is ri + 1. The rate at which these particles could move to the right is the sum of the persistence probability and the neighbor-weight probability, i.e., a+(1 a)( +i 1)=( + i + i 1). The other way for a particle to be located in bin i and move to the right is for a left-moving particle to choose to move to the right. This setup is shown in Fig. 4.4. Note that the total number of particles in each bin is the same and hence this setup produces the same neighbor weight as for a right-moving particle choosing to move to the right. These two options account for all possible ways for an extra particle to be present in bin i and move to the right to bin i + 1 with the system entering 77 the state f~r;~l; tg. Accordingly, the term ? in eq. (4.7) is given by ? : h a+ (1 a) +i 1 +i + i 1 i (ri + 1)P (J + i;r[~r];~l; t) + h (1 a) +i 1 +i + i 1 i (li + 1)P (J + i;l[~r;~l]; t): To explain ? in eq. (4.7), we need to account for particles that are present in bin i, move right to the bin i + 1, and leave the state f~r;~l; tg. In this setup, the neighbor-weight probability cri is exactly what we expect: (1 a)( + i )=( + i + i ). Again we attain the expressions in item ? by taking the product of the probability of the system being setup in the appropriate state, P (~r;~l; t), the number of particles available to move, either ri or li, and the probabilistic rate of particles moving to the right. This rate is the neighbor weight probability for both particles, with the persistence probability a added for right-moving particles. Hence, ? : h a+ (1 a) +i +i + i i riP (~r;~l; t) + h (1 a) +i +i + i i liP (~r;~l; t): By following similar reasoning, we can express fi in eq. (4.7), accounting for Right-moving particles : : : ri 1 ri ri+1 1 ri+2 : : : Left-moving particles : : : li 1 li + 1 li+1 li+2 : : : Total : : : ni 1 ni + 1 ni+1 1 ni+2 : : : Figure 4.4: Setup of particles for J+i;l[~r;~l] 78 particles entering state f~r;~l; tg when a particle in bin i moves to the left as fi : h a+ (1 a) i 1 +i + i 1 i (li + 1)P (~r; J i;l[~l]; t) + h (1 a) i 1 +i + i 1 i (ri + 1)P (J i;r[~r;~l]; t); and for fl in eq. (4.7), accounting for particles leaving the state f~r;~l; tg when a particle in bin i moves to the left, as fl : h a+ (1 a) i +i + i i liP (~r;~l; t) + h (1 a) i +i + i i riP (~r;~l; t): When the boundary conditions are periodic, the term in eq. (4.7) has an obvious structure. Adjustments should be made for di erent types of boundary conditions. 4.2.2 Deriving a system of ODEs Using the RDME (4.7), we de ne new variables which represent the average state of the system. We let Ri(t) = hri(t)i = P ~r P ~l riP (~r; ~l; t), where h i denotes expected value. Similarly, we de ne Li(t) = hli(t)i and Mi(t) = hni(t)i. By taking the partial derivative with respect to time t of each Ri and Li and by following the methodology of [5] and [29], we develop a system of ordinary di erential equations which describe our system. The resulting di erential equation for Ri, the right- moving particles in bin i, consists of particles entering from bin i 1 and exiting to 79 the right or to the left. Denoting the transition rates with T , the ODE for Ri is: dRi dt = hri 1T + ri 1i+ hli 1T + li 1 i hriT + ri i hriT ri i: (4.12) The equation for Li, the left-moving particles in bin i, consists of particles that enter from bin i+ 1 and exit to the left or right: dLi dt = hri+1T ri+1i+ hli+1T li+1 i hliT + li i hliT li i: (4.13) The transition rates are T+ri = c r i + a; (4.14) T+li = c r i ; (4.15) T li = c l i + a; (4.16) T ri = c l i; (4.17) where cri = (1 a) + i =( i + + i ) and c l i = (1 a) i =( i + + i ). Note that the rate at which particles move to the left or to the right depends on whether a particle is left-moving or right-moving. As de ned, the transition rate T+ri is the ?Rate of right-moving particles moving to the right? and the transition rate T li is the ?Rate of left-moving cells moving to the left?. Substituting eqs. (4.14){(4.17) 80 in (4.12){(4.13), we have dRi dt = aRi 1 Ri + hni 1c r i 1i; (4.18) dLi dt = aLi+1 Li + hni+1c l i+1i: (4.19) Observe that the di erential equations (4.18){(4.19) contain terms with an expected value. Our computational treatment of these terms will be discussed in Section 4.4. 4.3 An extended one-dimensional model of local interactions We now want to relax the assumption that all particles move at every time step. We still allow particles to persist in their direction of motion with probability a, and particles use neighbor-based weight to determine the probability of choosing the direction of motion; however, we now allow a particle or remain stationary with probability b. Hence, for a particle i with position xi(t) at time t, we calculate the new position at the next time step according to xi(t+ 1) = 8 >>>>>>>>>>>>>>< >>>>>>>>>>>>>>: xi(t) + pi(t); with probability a; xi(t); with probability b; xi(t) + x; with probability +i +i + i (1 a b); xi(t) x; with probability i +i + i (1 a b); (4.20) 81 where pi(t) = xi(t k) xi(t k 1) is the previous direction in which the particle moved with k = min[argmaxjjxi(t j) xi(t j 1)j], the length of time since the particle last changed its bin position. 4.3.1 Derivation of the Reaction-Di usion Master Equation The number of particles ni(t) in bin i at time t can be calculated from our previous system of N particles as ni(t) = NX j=1 xj(t);i: (4.21) At every time step, each cell has a preferred direction. In this extended model, an individual particle may stop. Consequently, there are four types of particles: right-moving, left-moving, right-stationary, and left-stationary. A particle is con- sidered to be right- (or left- )stationary if it moved to the right (or to the left) to enter the current bin but remained stationary during the previous time step. Note that a stationary particle remembers its previous direction. This memory allows for continued persistence in the preferred direction when the particle resumes moving. Let ri(t) be the number of particles that moved into bin i from bin i 1 during the previous time step. Let rsi (t) denote the number of particles at time t which have become stationary, but previously moved to the right from bin i 1 to i. Similarly, we de ne li(t) as the left-moving particles in bin i at time t, which moved left from bin i+ 1 into bin i during the previous time step. Finally, we de ne lsi (t) to be the number of particles in bin i at time t which moved left in their most recent transition 82 between bins, but remained stationary during the previous time step. We assume that particles become (or remain) stationary with probability b. The governing set of rules are shown in Figure 4.5. ri ri s l i li s l i!1 ri+1 ci l ci l a + ci l a + ci l b b a + ci r a + ci r ci r ci r b b Figure 4.5: A depiction of the options of motion for particles out of a right-moving ri or left-moving li position, with associated probabilities for a system of particles which are either right-moving or left-moving and have the ability to become stationary with memory of a preferred direction. Not shown are the sources of right-moving and left-moving particles. The source of particles in ri is all particle types in bin i 1, and the source of particles in li is all particle types in bin i + 1. The associated probabilities of these sources can be extracted from the diagram. The probability a denotes persistence, the probability b denotes becoming (or staying) stationary and the probabilities cri and c l i denote choosing a direction based on neighbor weights. In order to write a reaction-di usion master equation describing how these populations evolve probabilistically in time, we consider the probability density function P (~r;~l; ~rs; ~ls; t) describing the likelihood of a system to be in a given state f~r;~l; ~rs; ~ls; tg. For transitions between these four populations, there is the additional possibility of particles becoming stationary or moving after having been stationary, making the RDME slightly more complicated. Accordingly, the equation is written 83 as @P @t (~r;~l; t) = k 1X i=1 h Probability that a particle moves out of bin i to the right to enter the state f~r;~l; ~rs; ~ls; tg | {z } ? Probability that a particle moves out of bin i to the right to leave the state f~r;~l; ~rs; ~ls; tg | {z } ? i + kX i=2 h Probability that a particle moves out of bin i to the left to enter the state f~r;~l; ~rs; ~ls; tg | {z } fi Probability that a particle moves out of bin i to the left to leave the state f~r;~l; ~rs; ~ls; tg | {z } fl i + kX i=1 Probability that a particle in bin i becomes stationary to enter or leave the state f~r;~l; ~rs; ~ls; tg | {z } + Boundary terms | {z } ? (4.22) To address each of the terms in (4.22), we need to de ne additional particle- transposing operators. We de ne these operators so that they each act on the state of the system, f~r;~l; ~rs; ~lsg. In Table 4.1, we de ne all operators, but only show the populations which are directly a ected by the operator. All other populations remain whatever they are in the state f~r;~l; ~rs; ~lsg, as this table accounts for all possible translations which are one cell movement away from the state of interest. Given the operators from Table 4.1, we can de ne the terms in the RDME (4.22). For ? in (4.22), we consider all the terms where particles move to the right with the system entering state f~r;~l; ~rs; ~ls; tg. Note that particles from each 84 Operator Result J+i;r : : : ; ri + 1; ri+1 1; : : : J+i;rs : : : ; ri+1 1; : : : ; rsi + 1; : : : J+i;l : : : ; ri+1 1; : : : ; li + 1; : : : J+i;ls : : : ; ri+1 1; : : : ; lsi + 1; : : : J i;r : : : ; ri + 1; : : : ; li 1 1; : : : J i;rs : : : ; li 1 1; : : : ; rsi + 1; : : : J i;l : : : ; li 1 1; li + 1; : : : J i;ls : : : ; li 1 1; : : : ; lsi + 1; : : : Jsi;r : : : ; ri + 1; : : : ; rsi 1; : : : Jsi;l : : : ; li + 1; : : : ; lsi 1; : : : Table 4.1: De nitions of particle-transposition operators acting on f~r;~l; ~rs; ~lsg. The rst subscript indicates the bin on which the operator is acting. The second subscript indicates the population type. The superscript indicates the direction in which that particle is moving, right (+), left (-), or (s) if the particle is becoming stationary. population-type in bin i are capable of moving right, based on the neighbor-weight probability, which is calculated based on the preceding state of the system. Further- more, right-moving particles ri and right-stationary particles rsi both persist in their preferred direction with probability a. Consequently, the term ? in (4.22) takes the form ? : h a+ (1 a) +i 1 +i + i 1 i (ri + 1)P (J + i;r[~r;~l; ~rs; ~ls]; t) + h (1 a) +i 1 +i + i 1 i (li + 1)P (J + i;l[~r;~l; ~r s; ~ls]; t) + h a+ (1 a) +i 1 +i + i 1 i (rsi + 1)P (J + i;rs [~r;~l; ~rs; ~ls]; t) + h (1 a) +i 1 +i + i 1 i (lsi + 1)P (J + i;ls [~r;~l; ~r s; ~ls]; t): For item ? in (4.22), we consider all terms where particles move right, with 85 the system leaving state f~r;~l; ~rs; ~ls; tg. Hence ? : h a+ (1 a) +i +i + i i riP (~r;~l; ~rs; ~ls; t) + h (1 a) +i +i + i i liP (~r;~l; ~rs; ~ls; t) + h a+ (1 a) +i +i + i i rsiP (~r;~l; ~rs; ~ls; t) + h (1 a) +i +i + i i lsiP (~r;~l; ~rs; ~ls; t): For item fi in (4.22), we consider all terms where particles move left, with the system entering the state f~r;~l; ~rs; ~ls; tg: fi : h (1 a) i 1 +i + i 1 i (ri + 1)P (J i;r[~r;~l; ~rs; ~ls]; t) + h a+ (1 a) i 1 +i + i 1 i (li + 1)P (J i;l[~r;~l; ~r s; ~ls]; t) + h (1 a) i 1 +i + i 1 i (rsi + 1)P (J i;rs [~r;~l; ~rs; ~ls]; t) + h a+ (1 a) i 1 +i + i 1 i (lsi + 1)P (J i;ls [~r;~l; ~r s; ~ls]; t): In item fl in (4.22), we consider all terms where particles move left, with the system leaving the state f~r;~l; ~rs; ~ls; tg: fl : h (1 a) i +i + i i riP (~r;~l; ~rs; ~ls; t) + h a+ (1 a) i +i + i i liP (~r;~l; ~rs; ~ls; t) + h (1 a) i +i + i i rsiP (~r;~l; ~rs; ~ls; t) + h a+ (1 a) i +i + i i lsiP (~r;~l; ~rs; ~ls; t): In item in (4.22), we account for the probability of moving particles be- coming stationary, both to enter and exit the state f~r;~l; ~rs; ~ls; tg. Particles become 86 stationary with probabilistic rate b. Hence : b(ri + 1)P (J s i;r[~r;~l; ~rs; ~ls]; t) + b(li + 1)P (J s i;l[~r;~l; ~rs; ~ls]; t) b h ri + li i P (~r;~l; ~rs; ~ls; t): Finally, the boundary terms ? in eq. (4.22) take an obvious form in the case of periodic boundary conditions. Proper adjustments should be made for other types of boundary conditions. 4.3.2 Deriving a system of ODEs We de ne the expected value of each population, e.g. Ri = hrii = X ~r X ~l X ~rs X ~ls riP (~r;~l; ~rs; ~ls; t): (4.23) Using this de nition, taking the derivative with respect to time, substituting the RDME as appropriate and reindexing many of the summations, we obtain a new system of ODEs that models the average behavior of the system. We denote the transition rates with T . The average behavior of the right-moving particles in bin i, Ri, consists of particles entering from bin i 1 and leaving either to the right, to the left, or to the stationary compartment. Consequently: dRi dt = hri 1T + ri 1i+ hli 1T + li 1 i+ hrsi 1T + ri 1i+ hl s i 1T + li 1 i hriT + ri i hriT ri i hriT s rii: (4.24) 87 Similarly, the average behavior of the left-moving particles in slot i, Li, consists of particles entering from bin i+ 1 and leaving either to the right, to the left, or to the stationary compartment: dLi dt = hri+1T ri+1i+ hli+1T li+1 i+ hrsi+1T ri+1i+ hl s i+1T li+1 i hliT + li i hliT li i hliT s lii: (4.25) The right-stationary population in bin i, Rsi , increases as right-moving particles becoming stationary and decreases as the stationary particles leave to the left or the right. dRsi dt = hriT s rii hr s iT + ri i hr s iT ri i: (4.26) The left-stationary population in bin i, Lsi , increases as left-moving particles becoming stationary and decreases as the stationary particles leave to the left or the right. dLsi dt = hliT s lii hl s iT + li i hlsiT li i: (4.27) The transition rates T+ri , T + li , T li , and T ri are rede ned from (4.14){(4.17) by substituting 1 a b for 1 a in each expression. We de ne the transition rate for particles going from a moving compartment to a stationary compartment as: T sri = T s li = b: (4.28) Using the expressions for transition rates, we can simplify some of the expected value operators. The population of right-moving particles Ri simpli es to a system 88 with right-moving and stationary particles in bin i 1 persisting into bin i and all particles in bin i 1 moving to the right with a neighbor-weighted probability cri 1. Additionally, all right-moving particles leave the system at every time step, either by persisting, becoming stationary, or choosing to move toward a neighboring bin: dRi dt = a(Ri 1 +R s i 1) + hni 1c r i 1i Ri: (4.29) The population of right-stationary particles Rsi consists of particles in Ri becoming stationary with stopping rate b and leaving the stationary population with rate 1 b: dRsi dt = bRi (1 b)R s i : (4.30) The population of left-moving particles Li simpli es to a system with left-moving and stationary particles in bin i+1 persisting into bin i with rate a and all particles in bin i+ 1 moving to the left with a neighbor-weighted probability cli+1. Additionally, all left-moving particles leave the system at every time step, either by persisting, by becoming stationary, or by choosing to move toward a neighboring slot: dLi dt = a(Li+1 + L s i+1) + hni+1c l i+1i Li: (4.31) The population of left-stationary particles Lsi consists of particles in Li becoming stationary with stopping rate b and leaving the stationary population with rate 1 b. dLsi dt = bLi (1 b)L s i (4.32) 89 Note that this system, (4.29){(4.32), conserves the total number of particles in the system. We can see this by summing the di erential equations (4.29){(4.32) over all bins and population-types. In our simulations, in the case that a particle has no neighbors, we set the neighbor-weight to zero to avoid divide by zero errors. This does not a ect the conservative nature of the system. If particles are not allowed to enter a slot with the neighbor-weight probability due to the absence of neighbors, then there are no particles in that bin. Hence, no particles can leave that bin either (at all time points, all moving-type particles must enter a new particle- type/position). In this way, setting the neighbor-weight to zero is the appropriate way to maintain the conservation of particles. This is typically not relevant at high particle densities but can be important as aggregations form, leaving sections of empty bins. 4.4 Simulations The solution of the system of ordinary di erential equations (4.29){(4.32) is simulated with various initial and boundary conditions. To start, we only con- sider periodic boundary conditions. Furthermore, the four-population system can be reduced to the two-population system by setting b = 0 and setting the initial stationary populations to zero. For this reason, we primarily consider the more complex system with four populations in simulations. Observe that the ordinary di erential equations (4.29){(4.32) contain unre- solved expected value terms. To deal with these terms in our simulations, we 90 approximate each expected value term by the sum, product and quotient of the expected value of each term contained therein; that is, we drop the expected value operator and replace each term with its expected value. For instance, hriri 1i would be estimated by RiRi 1. We do this in order to expedite computation time and eliminate the need to calculate the probability density function P (~r;~l; ~rs; ~ls; t) for all possible particle number and time combinations. In the upcoming discussion of simulations, we present a large number of im- ages, in most of which we do not display numbers along certain axes to aid in visual discernment. We instead o er one example which illustrates these axes with more details. In Figure 4.6, we show a sample simulation of model (4.29){(4.32), with persistence probability a = 0:5, stopping probability b = 0:1, and neighbor detec- tion distance D = 7. In Figure 4.6(a), the i-axis varies from 0 (at the left) to 100 (at the right) and the t-axis ranges from 0 (at the front) to 1000 (at the back). In Figure 4.6(b), the time is xed at 1000 and the i-axis ranges from 0 (at the left) to 100 (at the right). The vertical axis in each case is for Mi(t), the expected total number of particles in bin i at time t. 4.4.1 Comparing the ODEs to the stochastic particle model Our rst numerical simulation compares the ODEs with the stochastic particle model. Our results are shown in Figures 4.7 and 4.8. In both gures, we demonstrate three di erent observed patterns in order to show how the ODE model (4.29){(4.32) can replicate the patterns obtained with the stochastic particle system (4.20). The 91 (a) 0 25 50 75 100 0 200 400 600 800 1000 0 10 20 30 40 50 i a = 0.5, D = 7.0 t M i(t) (b) 0 25 50 75 100 0 5 10 15 20 25 30 35 40 45 50 a = 0.5, D = 7.0 i M i(1000 ) Figure 4.6: A sample simulation of model (4.29){(4.32) with model parameters a = 0:5, b = 0:1 and D = 7. In (a), we show the evolution of the system in 100 bins over 1000 time steps. In (b), we show the state of the system at t = 1000. three patterns shown are the formation of aggregations of particles, the lack of aggregations, and the formation of aggregations that merge. In Fig. 4.7, we show three-dimensional surface plots illustrating the number of particles in 100 bins from time 0 to 1000. In Fig. 4.8, we o er the same simulations as in Fig. 4.7, except that the number of particles in bin is only portrayed at the nal time 1000. This allows us to observe the evolution of each system as well as the nal state of each simulation after 1000 time steps. Note that the agent based model (4.20) does not exactly replicate the predic- tions of the ODEs. To match the behavior of the ODEs, one must perform multiple stochastic simulations and then average over the results. In comparing the aggrega- tion patterns produced by the particle system to the ODE, we see a similar number of aggregations forming at approximately the same time. In the rst example with persistence parameter a = 0:3 and neighbor detection distance D = 10, four ag- 92 Aggregations i t 0 50 100 Deterministic model a = 0.3, D = 10.0 i t 0 100 200 Deterministic model a = 0.3, D = 20.0 No distinct aggregations i t 0 10 20 Deterministic model a = 0.9, D = 2.0 i t 0 10 20 Deterministic model a = 0.7, D = 7.0 Merging aggregations i t 0 50 Deterministic model a = 0.5, D = 7.0 i t 0 50 100 Deterministic model a = 0.5, D = 10.0 Figure 4.7: A comparison of the stochastic particle system to the deterministic ODE for di erent scenarios: forming aggregations, the absence of aggregations, and merging aggregations. The parameters a, the persistence probability, D, the neighbor detection distance, and the initial data are identical in each case. The parameters a and D are varied as indicated above each plot to capture the di erent scenarios. For each of these images, the stopping probability b is 0.1 and the number of bins k is 100. On these 3-D images, the i-axis is for bin number, from 1 to 100, the t-axis is for time, from 1 to 1000, and the vertical z-axis is for particle number. 93 Aggregations i 0 20 40 60 80 Deterministic model a = 0.3, D = 10.0 i 0 20 40 60 80 Stochastic model a = 0.3, D = 10.0 i 0 50 100 150 Deterministic model a = 0.3, D = 20.0 i 0 50 100 150 Stochastic model a = 0.3, D = 20.0 No distinct aggregations i 8 8.2 8.4 8.6 Deterministic model a = 0.9, D = 2.0 i 0 5 10 15 20 Stochastic model a = 0.9, D = 2.0 i 7.5 8 8.5 9 Deterministic model a = 0.7, D = 7.0 i 0 5 10 15 20 Stochastic model a = 0.7, D = 7.0 Merging aggregations i 0 20 40 60 Deterministic model a = 0.5, D = 7.0 i 0 20 40 60 80 Stochastic model a = 0.5, D = 7.0 i 0 20 40 60 80 Deterministic model a = 0.5, D = 10.0 i 0 20 40 60 Stochastic model a = 0.5, D = 10.0 Figure 4.8: Comparison of the stochastic particle system and the deterministic ODE for di erent scenarios: aggregations, the absence of aggregations, and merging aggregations at t = 1000. The parameters a, the persistence probability, D, the neighbor detection distance, and the initial data are the same in each case. The parameters a and D are varied as indicated above each plot to create the di erent scenarios. In all images, the stopping probability b is 0.1 and the number of slots k is 100. On these 2-D plots, the i-axis is for bin number, ranging from 1 to 100 and the vertical z-axis is for the number of particles, speci cally Mi(1000). 94 gregations are formed over the 100 bins in both the stochastic and deterministic simulation. In Fig. 4.8, we see that the forming aggregations do not end in identical locations and they are not the same height. Still, the patterns are similar: their number is similar and they appear to be of approximately the same width. Similarly, for the parameter set with persistence probability a = 0:3 and neighbor detection distance D = 20, we observe that the aggregation patterns match in number, width, and approximate height and location. In the case where no aggregations form, particles are seen to be distributed relatively uniformly, as observed in the 3-D surface evolution plots in Fig. 4.7. Upon closer inspection of this distribution, as can be observed in Fig. 4.8, the particles actually form a slight wave, the shape of which depends on the parameters and the initial conditions. We initially expected that simulations without aggregations would produce a uniform distribution of particles; however, this only appears to be the case when a + b = 1 and the initial data is uniformly distributed. We discuss this in more detail in the next section. In comparing the deterministic model and the stochastic model, we see that the stochastic model contains a lot of variation. The overall wave trend of the stochastic simulation seems to match the deterministic simulation but only loosely and not in magnitude. These wild variations in particle number for each bin illustrate the relatively large stochastic e ect of the particle system. These e ects become more important when looking at simulations with unstable dynamics. In the last case of simulations with merging aggregations, the merging times and patterns are very sensitive to the exact number of particles in each location. 95 In comparing the surface evolutions in Fig. 4.7 of the stochastic simulations to the deterministic simulations, we see that the aggregations do not merge at the same time and the larger (or smaller) aggregations are not in the same locations. However, we do observe that the general emerging patterns of particles are somewhat similar at the nal time t = 1000, especially in the rst merging aggregations example with a = 0:5 and D = 7. In the second merging aggregation example with a = 0:5 and D = 10, it appears as though the stochastic simulation will eventually form two large peaks, with the smaller peaks merging. We further explore patterns of merging aggregations in Figure 4.9, where we show four di erent instances of a simulation of the stochastic particle system with the same parameters (a = 0:5; b = 0:1; D = 7) and initial conditions. In the simulation of the deterministic ODE model, there are initially four aggregations. Eventually, after approximately 900 time steps, two of these aggregations merge. Of the resulting three aggregations, there are two relatively tall peaks and one shorter peak. In constrast, in the four stochastic simulations, three simulations yield similar results with three peaks, two tall and one short, at t = 1000. The simulation with four peaks at t = 1000 appears as though aggregations would be likely to merge shortly after this time. In the stochastic model, the initial data seem to converge quickly to ve aggregations, two of which in each image merge almost immediately. Immediate merging of two small aggregations may explain how the taller peak in the deterministic model developed. While a wide variation in the exact positions and heights of aggregations is observed in many simulations, we also observe qualitatively similar general patterns, 96 i t 0 50 Deterministic model a = 0.5, D = 7.0 i 0 20 40 60 Deterministic model a = 0.5, D = 7.0 i 0 20 40 60 80 Stochastic model a = 0.5, D = 7.0 i 0 20 40 60 Stochastic model a = 0.5, D = 7.0 i 0 20 40 60 80 Stochastic model a = 0.5, D = 7.0 i 0 20 40 60 Stochastic model a = 0.5, D = 7.0 Figure 4.9: Four simulations of the stochastic model for the same parameter set, alongside the simulation produced by the deterministic ODE model. In all images, the persistence probability a is 0.5, the neighbor detection distance D is 7, the stopping probability b is 0.1, and the number of bins k is 100. On the 3-D plots, the i-axis is for bin number, ranging from 1 to 100, the t-axis is for time step, from 1 to 1000, and the vertical z-axis is for total particle number, Mi(t). On the 2-D plots, the i-axis is for bin number and the vertical z-axis is for total particle number, Mi(1000). such as aggregation number and peak width, emerging under similar conditions. Due to the wide variation in patterns and the length of time to run each stochastic simulation (varying from 2 to 10 minutes depending on our choice of parameters), we have chosen not to perform a Monte Carlo simulation to con rm that the average behavior of the stochastic simulation does in fact match the system of ODEs. We are able to extract critical information from the system without performing this exercise. 4.4.2 Uniform initial conditions Before considering more simulations of this model, let us make a few obser- vations about the 4-population model (4.29){(4.32). Suppose each bin contains the same number, x, of right-moving and left-moving particles and the same number, y, 97 of right-stationary and left-stationary particles. Upon substitution, we observe that the total number of particles, Mi(t) = Ri(t) +Rsi (t) +Li(t) +L s i (t), is not changing with respect to time. Furthermore, the moving and stationary populations are also at a steady state if x = 1 bb y. So we expect that for uniform initial conditions, with the same number of particles in each bin and equal left- and right-moving (and sta- tionary) particles, the total number of particles will remain constant with a uniform distribution. We check this steady state in the following numerical simulation. In Fig- ure 4.10, we illustrate simulations where each particle-type-bin combination con- tains 2 particles, for a total of 8 particles in each bin. The simulations are run for combinations of parameters varying the persistence probability a and stopping probability b through a set of values (0.1, 0.3, 0.5, 0.7, 0.9). Note that a + b 1 is a hard constraint on the system; that is, the particles can either persist, stop or choose a new direction, but they cannot perform more than one of these actions at a time. Observe that these simulations suggest that a uniform distribution of particles is an unstable steady state of the system. To demonstrate that the uniform steady state is unstable, we perform the calculations with tighter integration error tolerances. The simulations are per- formed using ODE45 in MATLAB. The default integration error tolerances, used in Fig. 4.10, are ?RelTol? = 1e-3 and ?AbsTol? = 1e-6. In Figure 4.11, the error tolerances have been adjusted by over two orders of magnitude; we use ?RelTol? = 1e-5 and ?AbsTol? = 1e-8. Comparing Fig. 4.10 and Fig. 4.11, we observe that the solution diverges more slowly for stricter error tolerances. More speci cally for 98 0.9 i t 7.9999 8.0000 8.0001 a = 0.1, b = 0.9 0.7 i t 7.5 8 8.5 a = 0.1, b = 0.7 i t 7.9999 8.0000 8.0001 a = 0.3, b = 0.7 0.5 i t 0 50 100 a = 0.1, b = 0.5 i t 7.5 8 8.5 a = 0.3, b = 0.5 i t 7.9999 8.0000 8.0001 a = 0.5, b = 0.5 0.3 i t 0 50 100 a = 0.1, b = 0.3 i t 0 50 a = 0.3, b = 0.3 i t 5 10 15 a = 0.5, b = 0.3 i t 7.9999 8.0000 8.0001 a = 0.7, b = 0.3 0.1 i t 0 50 100 a = 0.1, b = 0.1 i t 0 50 100 a = 0.3, b = 0.1 i t 0 50 a = 0.5, b = 0.1 i t 7.99 8 8.01 a = 0.7, b = 0.1 i t 7.9999 8.0000 8.0001 a = 0.9, b = 0.1 " b; a! 0.1 0.3 0.5 0.7 0.9 Figure 4.10: Exploration of parameter space for parameters a, the persistence prob- ability, and b, the stopping probability, for uniform initial conditions (Mi(0) = 8 for all i). Note that there is a hard constraint that a+ b 1. For cases where a+ b = 1, particles do not have the ability to choose to move in a new direction; they can only stop or persist in their initially preferred direction. In all cases, the interaction distance D is 10 and the number of bins k is 100. The i-axis is the bin number, ranging from 1 to 100, the t-axis is for time, from 1 to 1000, and the vertical z-axis is for total particle number, Mi(t). 99 0.9 i t 7.9999 8.0000 8.0001 a = 0.1, b = 0.9 0.7 i t 7.5 8 8.5 a = 0.1, b = 0.7 i t 7.9999 8.0000 8.0001 a = 0.3, b = 0.7 0.5 i t 0 50 100 a = 0.1, b = 0.5 i t 7.9 8 8.1 a = 0.3, b = 0.5 i t 7.9999 8.0000 8.0001 a = 0.5, b = 0.5 0.3 i t 0 50 100 a = 0.1, b = 0.3 i t 0 50 a = 0.3, b = 0.3 i t 7.9999 8.0000 8.0001 a = 0.5, b = 0.3 i t 7.9999 8.0000 8.0001 a = 0.7, b = 0.3 0.1 i t 0 50 100 a = 0.1, b = 0.1 i t 0 50 100 a = 0.3, b = 0.1 i t 0 50 a = 0.5, b = 0.1 i t 7.9999 8.0000 8.0001 a = 0.7, b = 0.1 i t 7.9999 8.0000 8.0001 a = 0.9, b = 0.1 " b; a! 0.1 0.3 0.5 0.7 0.9 Figure 4.11: Exploration of parameter space for parameters a, the persistence prob- ability, and b, the stopping probability, for uniform initial conditions (Mi(0) = 8 for all i). This gure has the same parameter setup and initial conditions as those in 4.10; however, the MATLAB integration tolerances are stricter in this Figure to illustrate the numerical instability of these simulations. Note that there is a hard constraint that a+b 1. For cases where a+b = 1, particles do not have the ability to choose to move in a new direction; they can only stop or persist in their initially preferred direction. In all cases, the interaction distance D is 10 and the number of bins k is 100. The i-axis is for bin number, ranging from 1 to 100, the t-axis is for time, from 1 to 1000, and the vertical z-axis is for particle number. 100 parameter sets (a; b) = (0.3, 0.3), (0.3, 0.5), (0.5, 0.1), and (0.5, 0.3), there is a no- ticeable delay in aggregation formation. Furthermore, for the parameter sets (a; b) = (0.1, 0.1) and (0.1, 0.3), we see a shift in the location of the aggregations. As expected, for parameter values where a+ b = 1, the resulting simulation still yields a constant uniform distribution. 4.4.3 Parameter analysis by simulation We consider the case of initial conditions that follow a Poisson distribution. We generated a set of initial data with a Poisson distribution with mean 2 for each population type and bin. We use the same set of initial data for each simulation. The total number of particles for each bin, as well as a break down of particle number by population, are shown in Fig. 4.12. The total number of particles in this example is 820. i 0 2 4 6 8 Initial conditions for right?moving particles i 0 2 4 6 8 Initial conditions for right?stationary particles i 0 2 4 6 8 Initial conditions for left?moving particles i 0 2 4 6 8 Initial conditions for left?stationary particles i 0 5 10 15 20 Initial conditions for all particle types Figure 4.12: Initial conditions for Figures 4.13{4.22. The rst four plots, for right- moving Ri(0), right-stationary Rsi (0), left-moving Li(0) and left-stationary L s i(0) are used to initiate the ODEs (and the stochastic systems in Fig. 4.7, 4.8 and 4.9. The last plot is the sum of these four plots, Mi(0) and is what is visible in subsequent images. Particle numbers are distributed with a Poisson distribution with mean 2 for each population-type in each bin i. The i-axis varies from 0 to 100. The total number of particles is 820. In comparing simulations for this initial data set, our objective is to develop a better understanding of the parameter space. That is, we hope to better understand 101 the e ect of each parameter and their interdependencies. We begin by considering the e ect of varying both the persistence probability a and the stopping probability b. Each probability can vary between 0 and 1. In varying both parameters, we are constrained by a + b 1: the probability of persisting, stopping, and changing directions based on neighbor-weights must add up to one. In Figure 4.13, we show how the spatial distribution of particles evolves over 1000 time steps. In Figure 4.14, we show the distribution of particles at the nal time 1000. In both gures, the parameters a and b are chosen from the set f0:1; 0:3; 0:5; 0:7; 0:9g. We observe that as the persistence probability a is increased, the number of aggregations decreases. A similar trend occurs with the stopping probability. As b increases, the number of aggregations decreases; however, this relationship appears to also depend on the value of a. The dependence of the number of aggregations on persistence probability a is stronger than the dependence on the stopping probability b. Note that for a = 0:1, for this set of parameters, the only e ect of increasing b appears to be in the increased width of the peaks, as shown in Fig. 4.14. The width of peaks also increases with increased persistence probability a. Note that for simulations with a + b = 1, we observe oscillating waves of cells in Fig. 4.13. We also observe, for images with a+ b = 1 in Fig. 4.14, that the magnitude of these waves decreases for an increase in persistence probability a. Further, note that a few of the simulations do not appear to have reached a steady state, e.g., parameter sets (a,b) = (0.3, 0.5) and (0.7, 0.1), where nonzero dips between peaks appear as a viable source of merging aggregations. We next study the relation between the parameters a, the persistence prob- 102 0.9 i t 0 10 20 a = 0.1, b = 0.9 0.7 i t 0 50 100 a = 0.1, b = 0.7 i t 0 10 20 a = 0.3, b = 0.7 0.5 i t 0 50 100 a = 0.1, b = 0.5 i t 0 50 a = 0.3, b = 0.5 i t 0 10 20 a = 0.5, b = 0.5 0.3 i t 0 50 100 a = 0.1, b = 0.3 i t 0 50 100 a = 0.3, b = 0.3 i t 0 20 40 a = 0.5, b = 0.3 i t 0 10 20 a = 0.7, b = 0.3 0.1 i t 0 50 100 a = 0.1, b = 0.1 i t 0 50 100 a = 0.3, b = 0.1 i t 0 50 100 a = 0.5, b = 0.1 i t 0 20 40 a = 0.7, b = 0.1 i t 0 10 20 a = 0.9, b = 0.1 " b; a! 0.1 0.3 0.5 0.7 0.9 Figure 4.13: Exploring the parameter space for parameters a, the persistence prob- ability, and b, the stopping probability. Note that a + b 1. For cases where a + b = 1, particles cannot change their direction; they can only stop or persist in their initially preferred direction. In all cases, the interaction distance D is 10 and the number of bins is 100. The i-axis is for the bin number, the t-axis is time, from 1 to 1000, and the vertical z-axis is for particle number. The initial conditions are given in Figure 4.12. 103 0.9 i 7 7.5 8 8.5 9 a = 0.1, b = 0.9 0.7 i 0 20 40 60 80 a = 0.1, b = 0.7 i 7.5 8 8.5 9 a = 0.3, b = 0.7 0.5 i 0 50 100 a = 0.1, b = 0.5 i 0 20 40 60 a = 0.3, b = 0.5 i 7.5 8 8.5 9 a = 0.5, b = 0.5 0.3 i 0 50 100 a = 0.1, b = 0.3 i 0 20 40 60 a = 0.3, b = 0.3 i 0 10 20 30 40 a = 0.5, b = 0.3 i 8 8.2 8.4 8.6 a = 0.7, b = 0.3 0.1 i 0 50 100 a = 0.1, b = 0.1 i 0 20 40 60 80 a = 0.3, b = 0.1 i 0 20 40 60 80 a = 0.5, b = 0.1 i 0 10 20 30 a = 0.7, b = 0.1 i 8 8.2 8.4 8.6 a = 0.9, b = 0.1 " b; a! 0.1 0.3 0.5 0.7 0.9 Figure 4.14: The nal distribution at time t = 1000 of the simulations shown in Figure 4.13. 104 ability, and D, the neighbor detection distance. To do this, we again allow a to take values in the set f0:1; 0:3; 0:5; 0:7; 0:9g. We consider two ranges of values for D: f2; 5; 10; 20; 40g and f1; 3; 5; 7; 9g. The rst range allows us to consider the large scale e ects of doubling the parameter D. The second range allows us to consider the e ect of varying parameter D by relatively small increments. The second range is more likely the biologically reasonable range to consider, but understanding the larger scale e ects of D is also important. We set the stopping probability b = 0:1. Results are shown in Figures 4.15{4.18. There are 1000 time steps in the surface evolution images. The snapshot of the distribution of cells is taken at the nal time t = 1000. The number of bins equals 100. Once again, increasing the persistence probability a clearly decreases the number of aggregations and increases the width of aggregation peaks. Increasing the neighbor detection distance D also decreases the number of aggregations. Interestingly, the neighbor detection distance D does not appear to a ect the width of aggregations. In Figs. 4.15 and 4.16, we notice that doubling the interaction distance yields approximately half as many aggrega- tions. In all gures (Figs. 4.15{4.18), when the persistence probability is a = 0:7, the system develops harmonic frequencies for a neighbor detection distance D that is between 1 and 5. When this occurs, the magnitude of the wave (in the absence of aggregations) becomes very small. Figures 4.17 and 4.18 also illustrate that the number of aggregations is very sensitive to D for very low values of D, but not very sensitive for values of D larger than 5. Finally, we explore the dependencies between the stopping probability b and the neighbor detection distance D. For these simulations, we set the persistence 105 40 i t 0 200 400 a = 0.1, D = 40.0 i t 0 200 400 a = 0.3, D = 40.0 i t 0 100 200 a = 0.5, D = 40.0 i t 0 50 100 a = 0.7, D = 40.0 i t 0 10 20 a = 0.9, D = 40.0 20 i t 0 100 200 a = 0.1, D = 20.0 i t 0 100 200 a = 0.3, D = 20.0 i t 0 50 100 a = 0.5, D = 20.0 i t 0 50 100 a = 0.7, D = 20.0 i t 0 10 20 a = 0.9, D = 20.0 10 i t 0 50 100 a = 0.1, D = 10.0 i t 0 50 100 a = 0.3, D = 10.0 i t 0 50 100 a = 0.5, D = 10.0 i t 0 20 40 a = 0.7, D = 10.0 i t 0 10 20 a = 0.9, D = 10.0 5 i t 0 50 100 a = 0.1, D = 5.0 i t 0 50 100 a = 0.3, D = 5.0 i t 0 50 a = 0.5, D = 5.0 i t 0 10 20 a = 0.7, D = 5.0 i t 0 10 20 a = 0.9, D = 5.0 2 i t 0 50 100 a = 0.1, D = 2.0 i t 0 50 a = 0.3, D = 2.0 i t 0 10 20 a = 0.5, D = 2.0 i t 0 10 20 a = 0.7, D = 2.0 i t 0 10 20 a = 0.9, D = 2.0 " D; a! 0.1 0.3 0.5 0.7 0.9 Figure 4.15: Exploring the parameter space for parameters a, the persistence prob- ability, and D, the neighbor detection distance. In all images, the stopping proba- bility b is 0.1 and the number of bins is 100. The i-axis is for the bin number, from 1 to 100, the t-axis is the time, from 1 to 1000, and the vertical z-axis is for the particle number. The initial conditions are given in Figure 4.12. 106 40 i 0 100 200 300 400 a = 0.1, D = 40.0 i 0 100 200 300 a = 0.3, D = 40.0 i 0 50 100 150 a = 0.5, D = 40.0 i 0 20 40 60 a = 0.7, D = 40.0 i 8 8.2 8.4 8.6 a = 0.9, D = 40.0 20 i 0 50 100 150 200 a = 0.1, D = 20.0 i 0 50 100 150 a = 0.3, D = 20.0 i 0 20 40 60 80 a = 0.5, D = 20.0 i 0 20 40 60 a = 0.7, D = 20.0 i 8 8.2 8.4 8.6 a = 0.9, D = 20.0 10 i 0 50 100 a = 0.1, D = 10.0 i 0 20 40 60 80 a = 0.3, D = 10.0 i 0 20 40 60 80 a = 0.5, D = 10.0 i 0 10 20 30 a = 0.7, D = 10.0 i 8 8.2 8.4 8.6 a = 0.9, D = 10.0 5 i 0 20 40 60 a = 0.1, D = 5.0 i 0 20 40 60 80 a = 0.3, D = 5.0 i 0 20 40 60 a = 0.5, D = 5.0 i 8.15 8.2 8.25 8.3 8.35 a = 0.7, D = 5.0 i 8 8.2 8.4 8.6 a = 0.9, D = 5.0 2 i 0 20 40 60 a = 0.1, D = 2.0 i 0 20 40 60 a = 0.3, D = 2.0 i 7.5 8 8.5 9 a = 0.5, D = 2.0 i 8.1996 8.1998 8.2 8.2002 8.2004 a = 0.7, D = 2.0 i 8 8.2 8.4 8.6 a = 0.9, D = 2.0 " D; a! 0.1 0.3 0.5 0.7 0.9 Figure 4.16: The nal distribution at time t = 1000 of the simulations shown in Figure 4.15. 107 9 i t 0 100 200 a = 0.1, D = 9.0 i t 0 50 100 a = 0.3, D = 9.0 i t 0 50 100 a = 0.5, D = 9.0 i t 0 10 20 a = 0.7, D = 9.0 i t 0 10 20 a = 0.9, D = 9.0 7 i t 0 50 100 a = 0.1, D = 7.0 i t 0 50 100 a = 0.3, D = 7.0 i t 0 50 a = 0.5, D = 7.0 i t 0 10 20 a = 0.7, D = 7.0 i t 0 10 20 a = 0.9, D = 7.0 5 i t 0 50 100 a = 0.1, D = 5.0 i t 0 50 100 a = 0.3, D = 5.0 i t 0 50 a = 0.5, D = 5.0 i t 0 10 20 a = 0.7, D = 5.0 i t 0 10 20 a = 0.9, D = 5.0 3 i t 0 50 100 a = 0.1, D = 3.0 i t 0 50 a = 0.3, D = 3.0 i t 0 20 40 a = 0.5, D = 3.0 i t 0 10 20 a = 0.7, D = 3.0 i t 0 10 20 a = 0.9, D = 3.0 1 i t 0 50 100 a = 0.1, D = 1.0 i t 0 10 20 a = 0.3, D = 1.0 i t 0 10 20 a = 0.5, D = 1.0 i t 0 10 20 a = 0.7, D = 1.0 i t 0 10 20 a = 0.9, D = 1.0 " D; a! 0.1 0.3 0.5 0.7 0.9 Figure 4.17: Exploring the parameter space for parameters a, the persistence proba- bility, andD, the neighbor detection distance. In all images, the stopping probability b is 0.1 and the number of bins is 100. The i-axis is for bin number, ranging from 1 to 100, the t-axis is the time, from 1 to 1000, and the vertical z-axis is for the particle number. The initial conditions are given in Figure 4.12. 108 9 i 0 50 100 150 a = 0.1, D = 9.0 i 0 20 40 60 80 a = 0.3, D = 9.0 i 0 20 40 60 a = 0.5, D = 9.0 i 0 5 10 15 20 a = 0.7, D = 9.0 i 8 8.2 8.4 8.6 a = 0.9, D = 9.0 7 i 0 50 100 a = 0.1, D = 7.0 i 0 20 40 60 80 a = 0.3, D = 7.0 i 0 20 40 60 a = 0.5, D = 7.0 i 7.5 8 8.5 9 a = 0.7, D = 7.0 i 8 8.2 8.4 8.6 a = 0.9, D = 7.0 5 i 0 20 40 60 a = 0.1, D = 5.0 i 0 20 40 60 80 a = 0.3, D = 5.0 i 0 20 40 60 a = 0.5, D = 5.0 i 8.15 8.2 8.25 8.3 8.35 a = 0.7, D = 5.0 i 8 8.2 8.4 8.6 a = 0.9, D = 5.0 3 i 0 20 40 60 a = 0.1, D = 3.0 i 0 20 40 60 a = 0.3, D = 3.0 i 0 10 20 30 a = 0.5, D = 3.0 i 8.199 8.2 8.201 8.202 8.203 a = 0.7, D = 3.0 i 8 8.2 8.4 8.6 a = 0.9, D = 3.0 1 i 0 20 40 60 a = 0.1, D = 1.0 i 7 7.5 8 8.5 9 a = 0.3, D = 1.0 i 8.15 8.2 8.25 a = 0.5, D = 1.0 i 8.196 8.198 8.2 8.202 8.204 a = 0.7, D = 1.0 i 8 8.2 8.4 8.6 a = 0.9, D = 1.0 " D; a! 0.1 0.3 0.5 0.7 0.9 Figure 4.18: The nal distribution at time t = 1000 of the simulations shown in Figure 4.17. 109 probability a to 0:3 which results in a more interesting dynamics compared with the case when a = 0:1. The stopping probability b takes values in f0:1; 0:3; 0:5; 0:7g. Two sets of values are considered for D: f2; 5; 10; 20; 40g and f1; 3; 5; 7; 9g. The results are shown in Figures 4.19{4.22. Clearly, for the larger range of D values, shown in Figs. 4.19{4.20, the stopping probability b only a ects the simulations when b 0:5. This is most likely due to our choice of the persistence probability parameter a. Once again it is observed that if the detection distance D is doubled, the number of aggregations is approximately halved. For b = 0:7, where a + b = 1, we observe the same wave pattern to what was seen in Figs. 4.13{4.18. Similar results are obtained for low values of D, as can be seen in Figs. 4.21 and 4.22. There are a few interesting artifacts of these simulations that may be related to our choice of initial data set. The aggregations appear in approximately the same location in every simulation. This is very notable for simulations with one or two peaks. Additionally, in cases of two aggregations, the relative size of the peaks almost always decreases with increasing bin number. Similarly, for cases with four aggregations, the relative height of the peaks increases with an increasing bin number. Furthermore, for cases with no aggregations, the resulting wave is in approximately the same phase for each simulation, although the magnitude varies. Comparing these situations to the set produced with uniform initial data, the latter yields shifted sets of peaks for di ering simulations. 110 40 i t 0 200 400 b = 0.1, D = 40.0 i t 0 100 200 b = 0.3, D = 40.0 i t 0 100 200 b = 0.5, D = 40.0 i t 0 10 20 b = 0.7, D = 40.0 20 i t 0 100 200 b = 0.1, D = 20.0 i t 0 50 100 b = 0.3, D = 20.0 i t 0 50 100 b = 0.5, D = 20.0 i t 0 10 20 b = 0.7, D = 20.0 10 i t 0 50 100 b = 0.1, D = 10.0 i t 0 50 100 b = 0.3, D = 10.0 i t 0 50 b = 0.5, D = 10.0 i t 0 10 20 b = 0.7, D = 10.0 5 i t 0 50 100 b = 0.1, D = 5.0 i t 0 50 100 b = 0.3, D = 5.0 i t 0 20 40 b = 0.5, D = 5.0 i t 0 10 20 b = 0.7, D = 5.0 2 i t 0 50 b = 0.1, D = 2.0 i t 0 20 40 b = 0.3, D = 2.0 i t 0 10 20 b = 0.5, D = 2.0 i t 0 10 20 b = 0.7, D = 2.0 " D; b! 0.1 0.3 0.5 0.7 Figure 4.19: Exploring the parameter space for parameters b, the persistence prob- ability, and D, the neighbor detection distance. In all images, the stopping proba- bility a is 0.3 and the number of bins is 100. The i-axis is for the bin number, from 1 to 100, the t-axis is for the time, from 1 to 1000, and the vertical z-axis is for the particle number. The initial conditions are given in Figure 4.12. 111 40 i 0 100 200 300 b = 0.1, D = 40.0 i 0 50 100 150 200 b = 0.3, D = 40.0 i 0 50 100 150 b = 0.5, D = 40.0 i 7.5 8 8.5 9 b = 0.7, D = 40.0 20 i 0 50 100 150 b = 0.1, D = 20.0 i 0 50 100 b = 0.3, D = 20.0 i 0 20 40 60 b = 0.5, D = 20.0 i 7.5 8 8.5 9 b = 0.7, D = 20.0 10 i 0 20 40 60 80 b = 0.1, D = 10.0 i 0 20 40 60 b = 0.3, D = 10.0 i 0 20 40 60 b = 0.5, D = 10.0 i 7.5 8 8.5 9 b = 0.7, D = 10.0 5 i 0 20 40 60 80 b = 0.1, D = 5.0 i 0 20 40 60 b = 0.3, D = 5.0 i 0 10 20 30 40 b = 0.5, D = 5.0 i 7.5 8 8.5 9 b = 0.7, D = 5.0 2 i 0 20 40 60 b = 0.1, D = 2.0 i 0 10 20 30 40 b = 0.3, D = 2.0 i 7.5 8 8.5 9 b = 0.5, D = 2.0 i 7.5 8 8.5 9 b = 0.7, D = 2.0 " D; b! 0.1 0.3 0.5 0.7 Figure 4.20: The nal distribution at time t = 1000 of the simulations shown in Figure 4.19. 112 9 i t 0 50 100 b = 0.1, D = 9.0 i t 0 50 100 b = 0.3, D = 9.0 i t 0 50 b = 0.5, D = 9.0 i t 0 10 20 b = 0.7, D = 9.0 7 i t 0 50 100 b = 0.1, D = 7.0 i t 0 50 100 b = 0.3, D = 7.0 i t 0 20 40 b = 0.5, D = 7.0 i t 0 10 20 b = 0.7, D = 7.0 5 i t 0 50 100 b = 0.1, D = 5.0 i t 0 50 100 b = 0.3, D = 5.0 i t 0 20 40 b = 0.5, D = 5.0 i t 0 10 20 b = 0.7, D = 5.0 3 i t 0 50 b = 0.1, D = 3.0 i t 0 50 b = 0.3, D = 3.0 i t 0 10 20 b = 0.5, D = 3.0 i t 0 10 20 b = 0.7, D = 3.0 1 i t 0 10 20 b = 0.1, D = 1.0 i t 0 10 20 b = 0.3, D = 1.0 i t 0 10 20 b = 0.5, D = 1.0 i t 0 10 20 b = 0.7, D = 1.0 " D; b! 0.1 0.3 0.5 0.7 Figure 4.21: Exploring the parameter space for parameters b, the persistence prob- ability, and D, the neighbor detection distance. In all images, the stopping proba- bility a is 0.3 and the number of slots k is 100. The i-axis is for bin number, from 1 to 100, the t-axis is for time, from 1 to 1000, and the vertical z-axis is for the particle number. The initial conditions are given in Figure 4.12. 113 9 i 0 20 40 60 80 b = 0.1, D = 9.0 i 0 20 40 60 b = 0.3, D = 9.0 i 0 20 40 60 b = 0.5, D = 9.0 i 7.5 8 8.5 9 b = 0.7, D = 9.0 7 i 0 20 40 60 80 b = 0.1, D = 7.0 i 0 20 40 60 80 b = 0.3, D = 7.0 i 0 10 20 30 40 b = 0.5, D = 7.0 i 7.5 8 8.5 9 b = 0.7, D = 7.0 5 i 0 20 40 60 80 b = 0.1, D = 5.0 i 0 20 40 60 b = 0.3, D = 5.0 i 0 10 20 30 40 b = 0.5, D = 5.0 i 7.5 8 8.5 9 b = 0.7, D = 5.0 3 i 0 20 40 60 b = 0.1, D = 3.0 i 0 20 40 60 b = 0.3, D = 3.0 i 7 7.5 8 8.5 9 b = 0.5, D = 3.0 i 7.5 8 8.5 9 b = 0.7, D = 3.0 1 i 7 7.5 8 8.5 9 b = 0.1, D = 1.0 i 7.5 8 8.5 9 b = 0.3, D = 1.0 i 8.1 8.15 8.2 8.25 8.3 b = 0.5, D = 1.0 i 7.5 8 8.5 9 b = 0.7, D = 1.0 " D; b! 0.1 0.3 0.5 0.7 Figure 4.22: The nal distribution at time t = 1000 of the simulations shown in Figure 4.21. 114 4.5 Discussion The focus of this chapter was the development of a one-dimensional ODE model from the original stochastic model, proposed in Chapter 2. After making mi- nor changes to the original model, we began by developing a simple two-population model which only accounted for right-moving and left-moving particles. We then allowed the particles to become stationary, adding two additional populations to the model: right-stationary and left-stationary. After deriving this four-population model, we were able to more e ectively and e ciently explore the parameter space and develop a better understanding of our system. We began by comparing simulations from our stochastic model to the ODE model (4.29){(4.32). These simulations illustrated that the general trends in both models are the same. These trends include the number of aggregations, the width and height of these peaks, and the dynamics of merging aggregations. For a constant uniform initial distribution of particles, we showed that while the ODE system is supposed to be at steady state, the numerical simulations sug- gested that this is an unstable steady state. This reveals itself in numerical sim- ulations where integration error accumulates, eventually causing the presence of aggregations, though we expect a constant uniform steady state. Changing the nu- merical integration error tolerances has a small e ect on the temporal appearance and location of aggregations produced by small disturbances due to integration error. We analyzed the parameter space by considering a set of random initial condi- tions with a Poisson distribution, with which we were able to examine the complex 115 interplay between the persistence probability a, the stopping probability b, and the neighbor detection distance D. In considering the parameters individually, we note that increasing the persistence parameter a decreases the number of aggregations and increases the width of peaks. The e ect of increasing the stopping probability b is a decrease in the number of aggregations and an increase in the width of aggre- gations. We observed that increasing a and b e ectively decreases the probability of a particle being able to change directions in order to move toward a neighbor. In this way, particles are unable to form as tight of aggregations and instead form larger, less numerous peaks. Increasing the neighbor detection distance D decreases the number of peaks but does not a ect the width of peaks; instead, it appears to a ect the height of such peaks. Note that changing D does not change the total probability with which a particle can choose to move in a new direction, and a and b remain xed. In this way, the width of peaks is not a ected. Furthermore, doubling D appears to halve the number of these aggregations. We expect that for neighbor detection distances exceeding half the total number of available bins, there will ei- ther be one or no aggregations. In such a system, all particles are able to sense all other particles. We also observed that when particles are only capable of persisting in their current direction or stopping, maintained by the constraint a + b = 1, the deter- ministic result appear periodic in nature, with the shape of the curve depending on general shaped of the initial data. These results also occur for parameter sets where a+ b is close to one and the neighbor detection distance D is relatively small. The biological implications of this study motivate future research directions of 116 this work. We expect that the neighbor detection distance, which is possibly related to the length of pili on the surface of cyanobacteria or the di usion length scale of signaling molecules, is actually xed, or increases with time and varies from particle to particle. The stopping and persistence probabilities may also vary from particle to particle, and the values most likely take on a wide range, depending on temporal characteristics, e.g., how many polysaccharides the particles have produced and how many are in the agarose. These characteristics can be incorporated in future work. Mathematically, this model lends itself to further analysis, by way of derivation of a partial di erential equation. It is not immediately clear how to apply approximation theory in order to do this, as the neighbor detection distance could be allowed to remain xed or go to zero, while other ratios are held xed. It is likely that deriving such a model would produce an integro-di erential equation, with integrals being used to count neighbors within xed distances of the particle. Furthermore, the model can be extended to two dimensions. 117 Chapter 5 B7-H1 and Cancer Immune Interactions 5.1 Introduction B7-H1 is a surface protein which has been found on carcinomas of the lung, ovary, and colon and in melanomas but not on most normal tissues [25]. In the case of tumor immunity, experiments have shown that when the B7-H1 surface protein, also referred to as PD-L1 and CD274, is present on a tumor, cytotoxic T cells (CTLs) are less e ective at inducing apoptosis in the cancer cells [25, 41]. It is believed that B7-H1 forms a blockade against CTLs by interacting with PD-1 on the surface of CTLs [4]. A better understanding of the mechanism and dynamics may allow medical researchers to develop a cancer treatment schedule speci cally targeting this molecular shield, allowing the immune system to more e ectively repress a tumor. In this chapter, we present a dynamical model for the induction of apoptosis in cancer cells by cytotoxic T cells. We consider apoptosis by two di erent mechanisms: Fas/FasL binding and perforin. The model, developed in Section 5.2, is t to percent lysis data for B7-H1 transfected and mock protein transfected cancer cells exposed to cytotoxic T cells for periods of 4 and 12 hours. A formula for calculating percent lysis from cell population data is derived in Section 5.3. The results of our model tting and an analysis of the parameter values are presented in Section 5.4. We were able to show how our model can be used to t percent lysis data as well as 118 capture desirable population trends. The results of this chapter were published in the Bulletin of Mathematical Biology [33]. 5.2 A Mathematical model for CTL and tumor interaction In this section, we consider a model of cytotoxic T cells interacting with the immune system in vitro. We also give parameter estimates for the model. 5.2.1 Model In the experimental model of [41], activated cytotoxic T cells are cultured with cancer cells; hence, for this model we consider possible interactions between the two populations. In de ning the model, let C(t) denote the \uncomplexed" cancer cells at time t, X(t) denote complexes between cancer cells and CTLs at time t, and T (t) denote the total number of cancer cells at time t, i.e., C(t)+X(t). Let P (t) denote perforin in solution and let E(t) denote the CTLs or e ector cell population at time t. First consider the kinetic equation for complexes of cancer cells and e ector cells. Assuming mass action kinetics, complexes form at a rate proportional to the number of cancer cells C(t) and e ector cells E(t), i.e. k1C(t)E(t). They dissociate at a rate k2X(t). This association and dissociation via Fas/FasL binding are depicted in Figure 5.1. Complexes are also removed from the system when the CTL induces apoptosis in the tumor cell via Fas/FasL binding; we assume this happens at rate 119 + E XC Legend: FasL Fas Figure 5.1: Cartoon depiction of an e ector cell, or CTL, interacting with a cancer cell via Fas/FasL binding to form complex X k3X(t). Hence, the dynamic equation for X(t) reads: dX dt = k1C(t)E(t) k2X(t) k3X(t) (5.1) For the kinetic equation for perforin, we will assume that perforin is produced at a rate proportional to the number of e ector cells E(t), but only when there are very few e ector cells in the system. Otherwise, we assume that the rate of perforin production saturates with respect to E(t). To achieve this, we use Michaelis- Menten kinetics. We also allow the presence of cancer cells to inhibit the production of perforin if enough cancer cells are present; we account for this by dividing the production term by km2+C(t). Perforin is removed from the system when it interacts with cancer cells; this happens at rate k4C(t)P (t). The resulting equation is: dP dt = kpE(t) (km1 + E(t))(km2 + C(t)) k4C(t)P (t) (5.2) Finally, we consider the cancer cells which are not in complexes, for which we assume a logistic tumor growth with rate constants k and k5. Cells are removed 120 from this population when a complex forms at a rate of k1C(t)E(t) but are added back if the complex dissociates, which happens at a rate of k2X(t). Cells which undergo apoptosis via the Fas/FasL binding are eliminated. The other mechanism by which cells are induced to apoptosis is modeled by assuming mass action kinetics for interaction between cancer cells and perforin which we consider at a rate of k4C(t)P (t). All terms are combined into: dC dt = kC(t) k5C(t) 2 k1C(t)E(t) + k2X(t) k4C(t)P (t) (5.3) During the rst 12 hours of culturing the e ector and cancer cells, we assume the total number of e ector cells does not change; that is, E(t) +X(t) is constant. Since CTLs form complexes with cancer cells, the number of uncomplexed e ector cells E(t) at any time is calculated as E(0) X(t). This assumption is experimentally justi ed in [41]. Recall that the total number of cancer cells in solution T (t) is C(t) + X(t). In order to calculate percent lysis, we must compare this cancer cell population to one which has been cultured in the absence of CTLs, a population which we denote by T (t). As with the previous cancer cell population, we assume logistic growth, giving us the following equation. dT dt = kT (t) k5(T (t))2 (5.4) Note that the kinetic coe cients are the same for T (t) and the logistic growth 121 terms of C(t). The initial conditions are chosen to correspond to the experimental data in [41]. In each experiment, the initial concentration of P815 cancer cells is 104 cells per 200 L well. The cancer cells in each culture have either been transfected with B7-H1 or transfected with a mock surface protein. CTLs were added to each culture so that ratio of CTLs to cancer cells is xed at various e ector cell to target cell ratios. We denote these xed ratios in boldface as E=T. At the end of the culture period, either 4 or 12 hours, the percent lysis of the tumor cells by the CTLs is then determined by chromium release assay. We begin xing our initial conditions by scaling our variables so that C(0) = 1. The initial e ector cell concentration is then determined by a scale factor E=T. To correspond to the experimental data, values of E=T for time t of 4 hours are 0.25, 0.5, 1, 2, 4, 8, and 16. Values of E=T for time t of 12 hours are 0.08, 0.15, 0.25, 0.6, 1.25, and 2.5. In this way, the initial e ector cell concentration E(0) is equal to E=T times C(0). We assume that no Fas-FasL complexes X are present in the system at time zero, that is X(0) = 0. Hence, we also assume T (0) = C(0) + X(0) = 1. We assume that the initial amount of perforin in the system is proportional to the number of e ector cells for low levels of e ector cells but is constant for large levels of e ector cells. To do this, we use a Michaelis-Menten model for the initial condition as can be seen in Table 5.1. To check the validity of our model, we use percent lysis data found in [41]. 122 Description Equation \Uncomplexed" cancer cells C(0) = 1 Complexes of CTLs and cancer cells X(0) = 0 All cancer cells in system T (0) = 1 CTLs, or E ector cells E(0) = (E=T)C(0) Perforin P (0) = k6E(0) km+E(0) Cancer cells in absence of CTLs T (0) = 1 Table 5.1: Initial conditions for model Percent lysis can be calculated from our model: % lysis at time t = 100 T (t) T (t) T (t) (5.5) Our derivation of this expression and our assumptions can be found in x5.3. 5.2.2 Parameters In our model, we assume that all of the parameters, except for three, are the same when considering the interactions of both B7-H1+ and mock-transfected cancer cells with cytotoxic T cells. The presence of B7-H1 is expected to reduce the e ectiveness of CTLs in inducing apoptosis in the cancer cells by means of forming a blockade; hence, we allow k3, the rate parameter for apoptosis induced by the Fas- FasL mechanism and k4, the rate parameter for apoptosis induced by perforin and granzymes, to di er depending on the presence or absence of B7-H1. Experimental studies [27] and [45] have shown that in the presence of B7-H1, the production of perforin is inhibited. Assuming that this has to do with the extent of the blockade and the number of cancer cells present, we also allow km2, the Michaelis-Menten 123 term for perforin production inhibition by the presence of cancer cells, to di er based on the presence or absence of B7-H1. The parameter values utilized in model simulations can be found in Table 5.2. The values were obtained by minimizing the sum of square errors between the 4 hour and 12 hour percent lysis data and percent lysis as calculated by our model also at the 4 and 12 hour marks. The minimization was performed using the MATLAB function lsqnonlin. In order to obtain parameter values which t both the B7-H1+ and B7-H1 data, the sum of square errors being minimized was the weighted sum of square errors over both data sets. The weight for each term in the sum is the square of the inverse of the experimental error at that point. Note that the parameter values for k3 and km2 di er as expected. With k4 being the parameter for Fas/Fas-L induced cytolysis, it has been experimentally predicted that the presence of B7-H1 cancer cells would inhibit this pathway thus making the value of k3 for B7-H1+ less than that of k3 for B7-H1 . Hirano et al. suggest that their results imply that the perforin/granzyme pathway of cytolysis is not inhibited by B7-H1 [41]. Two other recent papers, [27] and [45], have experimentally shown that perforin production is greatly decreased in the presence of B7-H1+ cancer cells. In combination, these references imply that the initial amount of perforin within e ector cells can be used to induce apoptosis in both B7-H1+ and B7-H1 cancer cells, but after the rst 4 hours or so, the e ector cells in the presence of B7-H1 will not be able to produce as much perforin. This is re ected in our values of km2. It is interesting that the value of k4 for B7-H1+ is slightly larger than that k4 for B7-H1 . We will examine this further in the next section. 124 Param. B7-H1+ B7-H1 Interpretation Units k 0.035 cancer growth rate hr 1 k1 0.0001 complex formation 2 10 2 L/(cell hr) k2 0.0001 complex dissociation hr 1 k3 0.0001 190 apoptosis by complex hr 1 k4 3.0 2.2 apoptosis by perforin k5 0.003 k=Cancer carrying capacity 2 10 2 L/(cell hr) k6 0.63 initial ratio of P/E kp 0.097 ratio of P/E km 1 initial M-M term for P & E 50 cells/ L km1 2.2 M-M term for P & E 50 cells/ L km2 80 0.1 M-M term for P & C 50 cells/ L Table 5.2: Parameter values for model simulation 5.3 Expression for percent lysis Percent lysis experiments are done by Chromium release assay; the goal of such an experiment is to provide a scale by which to compare the e ectiveness of di erent lytic agents or e ector cells in lysing target cells. In order to conduct an experiment to determine percent lysis, the target cells must be incubated with radioactive chromium. The cells take in some of the chromium and the rest is washed away before the experiment begins. After some period of time t, percent lysis is calculated with the following formula: % lysis = 100 (experimental 51Cr release) (spontaneous 51Cr release) (maximum 51Cr release) (spontaneous 51Cr release) (5.6) Each of the variables in the equation for percent lysis (5.6) is calculated at time t. The spontaneous 51Cr release at time t is determined by assessing how much chromium is in solution at time t where the target cancer cells have not been in 125 the presence of e ector cells. This is a measure of how much chromium has been spontaneously released by the cells after time t. The maximum 51Cr release is determined by exposing the target cancer cells, again in the absence of e ector cells, to a lytic agent and then measuring the amount of chromium in the system. The experimental 51Cr release is determined by measuring the amount of chromium in solution in which target cells have been exposed to e ector cells for a time duration t. Mathematical formulae to predict percent lysis have been developed, e.g. see [77], but these models do not take into account the full dynamics of the system. To relate this to our model, we need to make a couple of assumptions. We assume that the target cells have a uniform concentration of chromium, and let (t) denote the amount of chromium present per viable cell at time t. We also assume that (t) is identical for the two target cell populations we are considering, i.e., those in the absence of e ector cells, T (t), and those in the presence of e ector cells, T (t). Noting that we previously assumed that target cells have the same growth dynamics and associated growth parameters in our model, it seems reasonable to assume that (t) is reduced only by spontaneous release and mitosis. If a target cell dies, by apoptosis or otherwise, this should not a ect the amount of chromium in the remaining viable cells. In this sense, we are assuming that once chromium leaves the cell the amount that di uses back into the cell is negligible. Using our assumptions, we can derive a new equation for percent lysis. The maximum amount of chromium released minus the amount released spontaneously is equal to the amount of chromium that is present in viable cells, which are in the 126 absence of e ector cells, that is, (maximum 51Cr release) (spontaneous 51Cr release) = (t)T (t) (5.7) A similar expression can be determined for the amount of chromium which is present in viable cells that are in the presence of e ector cells. If the maximum release minus the spontaneous release tells us how much chromium could possibly be present in the cells at time t and we know how much has been removed from the cells by lysis at time t using the experimental release minus the spontaneous release, then we can calculate the amount of chromium present in the cells at time t by taking the di erence between how much chromium could possibly be in the cells and how much has been removed from the viable cell population as can be seen below. [(max. 51Cr release) (spont. rel.)] [(exp. rel.) (spont. rel.)] = (t)T (t) (5.8) Combining equations (5.6), (5.7), and (5.8), we have obtained an equation for percent lysis that we can now use with our model variables: % lysis at time t = 100 T (t) T (t) T (t) (5.9) In the next section, we use this equation to compare our simulation results to avail- able experimental data. 127 5.4 Results and Discussion In Figure 5.2, we present the percent lysis calculation results of our model in comparison to experimental data. There are two noteworthy trends in the data. First, note that the experimental data at 4 hours is not that di erent for B7-H1+ and B7-H1 ; however after 12 hours the expected percent lysis for B7-H1+ is ap- proximately half of the B7-H1 . Secondly, it is noteworthy that as E=T increases, the experimental percent lysis appears to saturate. That is, even if the e ector cells are at a very high ratio to the cancer cells, the CTLs are not su ciently capable of suppressing the tumor in these ratios. This may imply the existence of another mechanism of cancer cell avoidance of lysis by CTLs. It is important to note that with only three parameters di ering between the B7-H1+ and B7-H1 equation sets, our model captures both of the these desirable trends. 4 hours 12 hours 2 4 6 8 10 12 14 160 20 40 60 80 E/T % lysi s Model B7?H1? Model B7?H1+ Experimental B7?H1? Experimental B7?H1+ 0.5 1 1.5 2 2.5 0 20 40 60 80 E/T % lysi s Figure 5.2: Percent lysis experimental data t by model at 4 and 12 hours for both B7-H1+ and B7-H1 for various e ector to target cell ratios. The experimental data is from [41]. In Figure 5.3, we observe the cell population dynamics over the course of 128 12 hours. In in vivo experiments in [41] lasting 25 to 50 days, it was observed that a tumor cell population which was transfected with a mock protein, instead of B7-H1, would increase in number initially, but then gradually decrease as CTLs managed to lyse the population. In a twelve hour in vitro experiment in [41], the mock transfected cancer cells gradually decreased in number while B7-H1-positive cancer cells linearly increased in number. In looking at our population dynamics, we also see this behavior as the total cancer cell population for B7-H1 cells cultured with cytotoxic T cells, T , decreases with time. The B7-H1+ cancer cell population decreases initially due to the large amount of perforin in the system, but then increases gradually. Note that the amount of perforin in the system, representative of the CTLs ability to induce apoptosis via this mechanism, decreases to a negligible amount for CTLs in the presence of B7-H1 but does not go to zero for CTLs in the absence of B7-H1. The complex dynamics are only signi cant in the presence of B7- H1 and are shown on the corresponding graph. The negligible amount of complexes present in the B7-H1 system can be accounted for by comparing the rates at which complexes form, dissociate, and induce apoptosis via Fas/FasL binding. The rate parameter for Fas k3 is much larger than both k1 and k2, implying that almost every complex that forms is almost immediately assigned to inducing apoptosis via the Fas/FasL mechanism. To ensure that we actually have attained reasonable parameter values for the rate of perforin-induced apoptosis k4, we performed a Latin hypercube sampling ex- periment varying B7-H1+ and B7-H1 k4. The parameter values for B7-H1 positive and negative systems were allowed to vary between 1 and 4. The value plotted in 129 B7-H1 B7-H1+ 0 2 4 6 8 10 12 0 1 2 3 4 5 Time (hours) 10 4 cells per 200 uL wel l T* T 10P 0 2 4 6 8 10 12 0 1 2 3 4 5 Time (hours) 10 4 cells per 200 uL wel l T* T 1000X 10P Figure 5.3: Population dynamics for cancer cells, transfected with either B7-H1+ or a mock protein in the presence of e ector cells T and in the absence of e ector cells T . The concentration of the Fas-FasL complexes is also plotted but very small relative to T . For the initial conditions, we chose E=T = 1. Figure 5.4 is the sum of errors squared for the plotted parameter set divided by the sum of errors squared for the parameter set in table 5.2. The errors are the di erence in percent lysis between the experimental values and the predicted values and the summation is over both 4 and 12 hours at their respective prescribed E=T ratios. Values close to one are depicted by black dots, while larger relative errors are shown with lighter dots. As can be seen in Figure 5.4, the values for perforin-induced apoptosis rates do not necessarily have to be di erent in order to get a reasonable, albeit suboptimal, t to the percent lysis data. This supports the assumption in [41] that B7-H1 does not a ect or decrease the rate of cytolysis by perforin. It is noteworthy that our model inaccurately treats perforin as a molecule with some uniform concentration in solution; perforin in vivo would be localized to neighborhoods between CTLs and cancer cells. However, given that the solution is cultured at relatively high cell concentrations in vitro with only e ector and target 130 1 2 3 41 1.5 2 2.5 3 3.5 4 k4 for B7?H1 ? k 4 for B 7? H 1+ 1 1.1 1.2 1.3 1.4 1.5 Figure 5.4: Latin hypercube sampling analysis of parameter k4, apoptosis rate for perforin interacting with cancer cells, where all other parameters are xed at values given in table 5.2. The value being plotted is the weighted sum of square errors for the plotted parameter values divided by the weighted sum of square errors for the parameter values in the table. cells present, a uniform perforin concentration is a reasonable rst approximation to a global average of activity due to perforin and granzymes. Regardless of the appli- cability of this assumption, the model does provide means of considering apoptosis by two separate mechanisms. These mechanisms are not mathematically equivalent under quasi-steady state approximations, that is assuming dX=dt and dP=dt are both zero. Such considerations of utilizing both apoptosis mechanisms are a likely key to the development of a more accurate and more mechanistic model describing CTL and cancer cell interaction. Considering that perforin-induced apoptosis is the primary mechanism of apop- 131 tosis during the time series considered, our results seem to agree with experimental data that B7-H1 has little e ect on the rate of perforin-induced apoptosis but does inhibit perforin production. In fact, if the parameters for k4 of B7-H1-positive and B7-H1-negative, are equated, the curve ts are achieved with comparable accuracy as is displayed in Figure 5.2. It is important to note that our assumptions about the rate limiting mechanisms, perforin production and Fas-FasL complex formation, might not be accurate. For instance, it is possible that apoptosis by these means is not rate limited by binding mechanism but is instead limited by signal transduction. Indeed, De Pillis et al. conclude in [24] that the law of mass action is not a su cient way to model cytolysis as a whole. Regarding the parameter estimates, considering that the data we used is strictly in vitro and we are interested in extending our model to t in vivo data, a more comprehensive analysis of the parameter space might be necessary. The formula which we derive in Section 5.3 for calculating percent lysis from cell population data is simple, intuitive, and computationally cheap. While the formula may not be as precise as it could be, the development and validation of a more involved model may require experimental research. It is also noteworthy that accuracy associated with the use of this formula will depend on the assumed model dynamics of cancer cells in the absence of CTLs. In our case, we modeled cancer cell growth and death with a logistic growth model. This assumption could possibly be improved by using a cancer growth model which more accurately ts a growth curve for P815 cancer cells grown in vitro; however, this data was unavailable for the experiment. 132 5.5 Conclusions We have provided a model for CTL-induced apoptosis of cancer cells which is more mechanistic than those currently available while still being relatively elemen- tary. We were able to show how the model can be used to t percent lysis data for in vitro interactions between CTLs and cancer cells, transfected with either B7-H1 or a mock protein. Only three of the parameters were assumed to be di erent: the rate parameter for Fas-mediated apoptosis, the rate parameter for perforin-induced apoptosis and a Michaelis-Menten-like parameter for perforin production inhibition by cancer cells. The model?s projected cell population data behaved as desired with the mock-transfected cancer cell population decreasing in size and the B7-H1- transfected cancer cell population increasing in size. While the model may require more work to ensure realistic parameter values, the idea of incorporating both mech- anisms of apoptosis and a formula for calculating percent lysis were presented. For future work, in order to make inferences with the model at time inter- vals longer than twelve hours, we intend to relax the assumption that the CTL population is constant and incorporate T cell dynamics. It would also be useful to incorporate data speci cally involving perforin, Fas/FasL, and cell knockouts of each. Additionally, it might be possible to reduce the dimension of the parameter space by making quasi-steady state approximations for the concentrations of per- forin and CTL-cancer complexes. We would also like to extend our model to in vivo data, including the upregulation of B7-H1 on tumor cells and interaction with an antibody; however, dynamic data for B7-H1 upregulation is not currently available. 133 Chapter 6 General conclusions In the rst part of this dissertation, we presented mathematical models and simulations of phototaxis and local interactions between cyanobacteria. In the sec- ond part of this dissertation, we studied the dynamics of B7-H1 and cancer immune interactions. Phototaxis is the ability of certain microorganisms to migrate toward light. This motion is typically associated with ngering patterns similar to those that develop in unstable fronts. We focused on a common freshwater microorganism Synechocystis sp., a cyanobacterium that has been studied in a laboratory setting in order to understand the functionality of the cell and how the motion of individual cells is translated into emerging patterns on macroscopic scales. Experimentally, it is observed that a critical number of cells are necessary for the detection of the direction of light and for initiating a directional movement. While signi cant advances have been made in studying the ways by which the bacteria can detect the light and move towards it, it is not well understood how these unicellular organisms identify the direction of the light. It is also not understood how the motion of the individual cells and the cell-to-cell interactions are then translated into the observed complex patterns. Using the experimental data of Synechocystis sp. provided by the Devaki lab, 134 we were able to identify a collection of rules for the local movement of cells on small length- and time-scales. Synechocystis displays a quasi-random motion in which cells do not all necessarily move in the direction of light and instead form temporary aggregations of cells. These observations were then used to develop a series of stochastic mathematical models, which integrate the local rules of motion (on small scales), with the global, large-scale forcing due to the light source. We also proposed a stochastic model of local interactions in one dimension, with which we are able to derive a system of ordinary di erential equations. Simulations of these mathematical models provide interesting insights on the relative role of the local and global forcing. They also provide a tool for studying the critical threshold that determines whether local aggregation patterns form or not. The simulations of our models that incorporate global forcing by light, support the hypothesis that not all cells are capable of independently identifying the direction of light. Our mathematical modeling results are shown to be in agreement with the available experimental data. Finally, we presented a mathematical model of cancer immune interactions, as related to B7-H1, an anti-apoptotic receptor on cancer cells. B7-H1 is found in carcinomas of the lung, ovary, and colon and in melanomas. In mouse model and in vitro experiments, B7-H1-positive cancer cells have been shown to be immune resis- tant, while B7-H1-negative tumor cells are signi cantly more susceptible to being repressed by the immune system [25]. In order to develop a better understanding of this complex system, we derived a system of ordinary di erential equations, model- ing the two possible mechanisms of apoptosis of cancer cells by the immune system. 135 We isolated three parameters that depend on the presence or absence of B7-H1. A better understanding of this relationship may allow medical researchers to develop a cancer treatment speci cally targeting this molecular shield, allowing the immune system to more e ectively repress a tumor. 136 Bibliography [1] M. S. Alber, M. A. Kiskowski, J. A. Glazier, and Y. Jiang. On cellular automa- ton approaches to modeling biological cells. In Mathematical Systems Theory in Biology, Communication, and Finance, D. N. Arnold and F. Santosa, editors (IMA 134, Springer-Verlag, New York, 2002), 1{40. [2] B. Alberts, A. Johnson, J. Lewis, M. Ra , K. Roberts, and P. Walter. Molecular Biology of the Cell, 5th edition. Garland Science, NY, 2008. [3] I. Aoki. A simulation study on the schooling mechanism in sh. Bulletin of the Japanese Society of Scienti c Fisheries 48(8) (1982), 1081{1088. [4] T. Azuma, S. Yao, Z. Gefeng, A. S. Flies, S. J. Flies, and L. Chen. B7-H1 is a ubiquitous antiapoptotic receptor on cancer cells. Blood 111(7) (2008), 3635{3643. [5] R. E. Baker, C. A. Yates, and R. Erban. From microscopic to macroscopic descriptions of cell migration on growing domains. Bulletin of Mathematical Biology 72 (2010), 719{762. [6] N. Bellomo, A. Bellouquid, and E. DeAngelis. The modelling of the immune competition by generalized kinetic (Boltzmann) models: review and research perspectives. Mathematical and Computer Modelling 37 (2003), 65{86. [7] K. Berg, P. K. Selbo, L. Prasmickaite, T. E. Tjelle, S. Kjolsrud, G. H. Rodal, H. Anholt, K. Sandvig, J. Moan, G. Gaudernack, O. Fodstad, S. K. Rodal, and A. Hogset. Photochemical internalisation. A novel technology for site-speci c delivery of macromolecules into cytosol. Cancer Research 59 (1999), 1180{1183. [8] R. Berger, R. Rotem-Yehudar, G. Slama, S. Landes, A. Kneller, M. Leiba, M. Koren-Michowitz, A. Shimoni, and A. Nagler. Phase I safety and phar- macokinetic study of CT-011, a humanized antibody interacting with PD-1, in patients with advances hematologic malignancies. Clinical Cancer Research 14(10) (2008), 3044{3051. [9] D. Bhaya, A. Takahashi, P. Shahi, and A. R. Grossman. Novel motility mutants of Synechocystis strain PCC 6803 generated by in vitro transposon mutagenesis. Journal of Bacteriology 183 (2001), 6140{6143. [10] D. Bhaya, A. Takahashi, and A. R. Grossman. Light regulation of type IV pilus- dependent motility by chemosensor-like elements in Synechocystis PCC6803. Proceedings of the National Academy of Sciences 98(13) (2001), 7540{7545. [11] D. Bhaya. Light matters: phototaxis and signal transduction in unicellular cyanobacteria. Molecular Microbiology 53 (2004), 745{754. 137 [12] D. Bhaya, K. Naksugi, F. Fazeli, and M. S. Burriesci. Phototaxis and impaired motility in adenylyl cyclase and cyclase receptor protein mutants of Synechocys- tis sp. strain PCC 6803. Journal of Bacteriology 188(20) (2006), 7306{7310. [13] D. Bhaya, D. Levy, and T. Requeijo. Group dynamics of phototaxis: Inter- acting stochastic many-particle systems and their continuum limit. Hyperbolic Problems: Theory, Numerics, Applications (2008), 145{159. [14] U. Burkhart and D. P. Hader. Phototactic attraction in light trap experiments: A mathematical model. Journal of Mathematical Biology 10 (1980), 257{269. [15] M. Burriesci and D. Bhaya. Tracking phototactic responses and modeling motil- ity of Synechnocystis sp. strain PCC6803. Journal of Photochemistry and Pho- tobiology B: Biology 91 (2008), 77{86. [16] L. Chen. Co-inhibitory molecules of the B7-CD28 family in the control of T-cell immunity. Nature Reviews Immunology 4 (2004), 336{347. [17] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks. Collective memory and spatial sorting in animal groups. Journal of Theoretical Biology 218 (2002), 1{11. [18] F. Cucker and S. Smale. On the mathematics of emergence. Japanese Journal of Mathematics 2 (2007), 197{227. [19] F. Cucker and S. Smale. Emergent behavior in ocks. IEEE Transactions on Automatic Control 52 (2007), 852{862. [20] E. De Angelis and P. E. Jabin. Qualitative analysis of a mean eld model of tumor-immune system competition. Mathematical Models and Methods in Applied Sciences 13 (2003), 187{206. [21] E. De Angelis and P. E. Jabin. Mathematical models of therapeutical actions related to tumour and immune system competition. Mathematical Methods in the Applied Sciences 28(17) (2005), 2061{2083. [22] P. Degond and S. Motsch. Continuum limit of self-driven particles with orien- tation interaction. Mathematical Models and Methods in Applied Sciences 18 (2008), 1193{1215. [23] J. Delgado and J. S. Gonzalez-Garcia. Evaluation of spherical shapes swim- ming e ciency at low Reynoles number with applications to some biological problems. Physica D 168{169 (2002), 365{378. [24] L. G. de Pillis, A. E. Radunskaya, and C. L. Wiseman. A validated mathemati- cal model of cell-mediated immune response to tumor growth. Cancer Research 65(17) (2005), 7950{7958. 138 [25] H. Dong, S. E. Strome, D. R. Salomao, H. Tamura, F. Hiaro, D. B. Flies, P. C. Roche, J. Lu, and G. Zhu. Tumor-associated B7-H1 promotes T-cell apoptosis: a potential mechanism of immune evasion. Nature Medicine 8(8) (2002), 793{ 800. [26] K. Drescher, R. E. Goldstein, and I. Tuval. Fidelity of adaptive phototaxis. Proceedings of the National Academy of Sciences 107(25) (2010), 11171{76. [27] K. Ebelt, G. Babaryka, B. Frankenberger, C. G. Stief, W. Eisenmenger, T. Kirchner, D.J. Schendel, and E. Noessner. Prostate cancer lesions are sur- rounded by FOXP3+, PD-1+ and B7-H1+ lymphocyte clusters. European Journal of Cancer 45(9) (2009), 1664{1672. [28] F. Entschladen, M. Gunzer, C. M. Scheu ele, B. Niggemann, and K. S. Znker. T lymphocytes and neutrophil granulocytes di er in regulatory signaling and migratory dynamics with regard to spontaneous locomotion and chemotaxis. Cellular Immunology 199(2) (2000), 104{114. [29] R. Erban, S. J. Chapman, and P. Maini. A practical guide to stochas- tic simulations of reaction-di usion processes. 34 pages, available as http://arxiv.org/abs/0704.1908 (2007). [30] H. Fatehi, M. Meyer-Hermann, and M. T. Figge. Modelling cellular aggregation induced by chemotaxis and phototaxis. Mathematical Medicine and Biology 27 (2010), 373{384. [31] A. Galante, D. Levy, and C. Tomasetti. A mathematical model for microenvi- ronmental control of tumor growth. In IFMBE Proceedings 32, K. E. Herold, W. E. Bentley, and J. Vossoughi, editors (Springer, College Park, MD, 2010): 213-216. [32] A. Galante, S. Wisen, D. Bhaya, and D. Levy. Stochastic models and simu- lations of phototaxis. In Unifying Themes in Complex Systems Volume VIII: Proceedings of the Eighth International Conference on Complex Systems. H. Sayama, A. A. Minai, D. Braha, D. and Y. Bar-Yam, editors (New England Complex Systems Institute Series on Complexity: NECSI Knowledge Press, Boston, MA, 2011), 105{119. [33] A. Galante, K. Tamada, and Doron Levy. B7-H1 and a mathematical model for cytotoxic T cell and tumor cell interaction. Bulletin of Mathematical Biology 74(1) (2012), 91{102. [34] E. N. Golovchenko, L. G. Hanin, S. H. Kaufmann, K. V. Tyurin, and M. A. Khanin. Dynamics of granzyme B-induced apoptosis: mathematical modeling. Mathematical Biosciences 212 (2008), 54{68. [35] S. Y. Ha and E. Tadmor. From particle to kinetic and hydrodynamic descrip- tions of ocking. Kinetic and Related Models 1(3) (2008), 415{435. 139 [36] S. Y. Ha, K. Lee, and D. Levy. Emergence of time-asymptotic ocking in a stochastic cucker-smale system. Communications in Mathematical Sciences 7 (2009), 453{469. [37] D. P. Hader and M. Lebert. Photoorientation in photosynthetic agellates. In Chemotaxis, T. Jin and D. Hereld, editors (Methods in Molecular Biology 571, Humana Press, 2009), 51{65. [38] K. Hall. Directed evolution of light-activated drugs. Medical Hypotheses 53(6) (1999), 504{506. [39] W. Haupt. Phototaxis in plants. In International Review of Cytology 19, G. H. Bourne, editor (1966), 267{299. [40] T. Hillen and K. J. Painter. A users guide to PDE models for chemotaxis. Journal of Mathematical Biology 57 (2009), 183-217. [41] F. Hirano, K. Katsumi, H. Tamura, H. Dong, S. Wang, M. Ichikawa, C. Rietz, D.B. Flies, J.S. Lau, G. Zhu, K. Tamada, and L. Chen. Blockade of B7-H1 and PD-1 by monoclonal antibodies potentiates cancer therapeutic immunity. Cancer Research 65(3) (2005), 1089{1096. [42] D. Horstmann. From 1970 until present: the Keller-Segel model in chemotaxis and its consequences I. Jahresberichte DMV 105(3), (2003), 103{165. [43] A. Huth and C. Wissel. The simulation of movement of sh schools. Journal of Theoretical Biology 156(3) (1992), 365{385. [44] M. Ikeuchi and T. Ishizuka. Cyanobacteriochromes: a new superfamily of tetrapyrrole-binding photoreceptors in cyanobacteria. Photochemical & Pho- tobiological Sciences 7 (2008), 1159{1167. [45] H. Y. Jeong, Y. J. Lee, S. K. Seo, S. W. Lee, S. J. Park, J. N. Lee, H. S. Sohn, S. Yao, L. Chen, and I. Choi. Blocking of monocyte-associated B7-H1 (CD274) enhances HCV-speci c T cell immunity in chronic hepatitis C infection. Journal of Leukocyte Biology 83(3) (2008), 755{764. [46] D. Kaiser. Myxococcus|from single-cell polarity to complex multicellular pat- terns. Annual Review of Genetics 42 (2008), 109{130. [47] E. F. Keller and L. A. Segel. Model for chemotaxis. Journal of Theoretical Biology 30(2) (1971), 225{234. [48] D. Kirschner and J. C. Panetta. Modeling immunotherapy of the tumor-immune interaction. Journal of Mathematical Biology 37 (1998), 235{252. [49] M. Kucia, K. Jankowski, R. Reca, M. Wysoczynski, L. Bandura, D. J. Al- lendorf, J. Zhang, J. Ratajczak, M. Z. Ratajczak. CXCR4-SDF-1 signalling, locomotion, chemotaxis and adhesion. Journal of Molecular Histology 3 (2004), 233245. 140 [50] H. W. Kuhlmann, R. Braucker, and A. G. Schepers. Phototaxis in Porpos- toma notatum, a marine scuticociliate with a composed crystalline organelle. European Journal of Protistology 33 (1997), 295{304. [51] V. A. Kuznetsov, I. A. Makalkin, M. A. Taylor, and A. S. Perelson. Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis. Bulletin of Mathematical Biology 56 (2) (1994), 295{321. [52] R. Lai and T. L. Jackson. A mathematical model of receptor-mediated apop- tosis: dying to know why FasL is a trimer. Mathematical Biosciences and En- gineering 1(2) (2004), 325{338. [53] W. C. Lam, E. J. Delikatny, F. W. Orr, J. Wass, J. Varani, and P. A. Ward. The chemotactic response of tumor cells. A model for cancer metastasis. American Journal of Pathology 104(1) (1981), 69{76. [54] D. Levy and T. Requeijo. Modeling group dynamics of phototaxis: from particle systems to PDEs. Discrete and Continuous Dynamical Systems - Series B 9 (2008), 103{128. [55] D. Levy and T. Requeijo. Stochastic models for phototaxis. Bulletin of Math- ematical Biology 70 (2008), 1684{706. [56] D. Levy and S. Y. Ha. Particle, kinetic and uid models for phototaxis. Discrete and Continuous Dynamical Systems{Series B 12 (2009), 77{108. [57] Y. X. Li, R. Lukeman, and L. Edelstein-Keshet. Minimal mechanisms for school formation in self-propelled particles. Physica D 237 (2008), 699{720. [58] R. Lukeman, Y. X. Li, and L. Edelstein-Keshet. Inferring individual rules from collective behavior. Proceedings of the National Academy of Sciences 107 (2010), 12576{12580. [59] A. F. M. Maree, A. V. Pan lov, and P. Hogeweg. Phototaxis during the slug stage of Dictyostelium discoideum: a model study. Proceedings of the Royal Society B 266 (1999), 1351{1360. [60] A. Mogilner and G. Oster. Cell motility driven by actin polymerization. Bio- physical Journal 71(6) (1996), 3030{3045. [61] A. Nacev, S. H. Kim, J. Rodriguez-Canales, M. Tangrea, B. Shapiro, and M. R. Emmert-Buck. A dynamic magnetic shift method to increase nanoparti- cle concentration in cancer metastases: A feasibility study using simulations on autopsy specimens. International Journal of Nanomedicine 2011:6 (2011), 2907{2923. [62] M. P. Neilson, J. A. Mackenzie, S. D. Webb, and R. H. Insall. Modeling cell movement and chemotaxis using pseudopod-based feedback. SIAM Journal on Scienti c Computing 33(3) (2011), 1035{1057. 141 [63] D. V. Nicolau, J. P. Armitage, and P. K Maini. Directional persistence and the optimality of run-and-tumble chemotaxis. Computational Biology and Chem- istry 33 (2009), 269{274. [64] F. Niyonsaba, K. Iwabuchi, A. Someya, M. Hirata, H. Matsuda, H. Ogawa, and I. Nagaoka. A cathelicidin family of human antibacterial peptide LL-37 induces mast cell chemotaxis. Immunology 106 (2002), 20{26. [65] H. G. Othmer and T. Hillen. The di usion limit of transport equations II: Chemotaxis equations. SIAM Journal of Applied Mathematics 62(4) (2002), 1222{1250. [66] J. K. Parrish, S. V. Viscido, and D. Grunbaum. Self-organized sh schools: An examination of emergent properties. Biological Bulletin 202 (2002), 296{305. [67] C. S. Patlak. Random walk with persistence and external bias. Bulletin of Mathematical Biophysics 15 (1953), 311{338. [68] E. T. Roussos, J. S. Condeelis, and A. Patsialou. Chemotaxis in cancer. Nature Reviews Cancer 11 (2011), 573{587. [69] R. W. Ruddon. Cancer Biology, 4th edition. Oxford University Press, 2007. [70] J.H. Russell and T.J. Ley. Lymphocyte-mediated cytotoxicity. Annual Review of Immunology 20 (2002), 323{370. [71] C. Sawyer, J. Sturge, D. C. Bennett, M. J. O?Hare, W. E. Allen, J. Bain, G. E. Jones, and B. Vanhaesebroeck. Regulation of breast cancer cell chemotaxis by the phosphoinositide 3-kinase p110 . Cancer Research 63 (2003), 1667{1675. [72] I. F. Sbalzarini and P Koumoutsakos. Feature point tracking and trajectory analysis for video imaging in cell biology. Journal of Structural Biology 151(2) (2005), 182{195. [73] P. K. Selbo, G. Sivam, O. Fodstad, K. Sandvig, and K. Berg. In vivo documen- tation of photochemical internalisation, a novel approach to site speci c cancer therapy. International Journal of Cancer 92 (2001), 761{766. [74] B. Shapiro. Towards dynamic control of magnetic elds to focus magnetic car- riers to targets deep inside the body. Journal of Magnetism and Magnetic Ma- terials 321(10) (2009), 1594{1599. [75] J. D. Shields, M. E. Fleury, C. Yong, A. A. Tomei, G. J. Randolph, and M. A. Swartz. Autologous chemotaxis as a mechanism of tumor cell homing to lymphatics via interstitial ow and autocrine CCR7 signaling. Cancer cell 11(6) (2007), 526{538. [76] J. Teran, L. Fauci, and M. Shelley. Viscoelastic uid response can increase the speed and e ciency of a free swimmer. Physical Review Letters 104 (2010), 038101. 142 [77] R. M. Thorn and C. S. Henney. Kinetic analysis of target cell destruction by e ector T cells: I. Delineation of parameters related to the frequency and lytic e ciency of killer cells. Journal of Immunology 117 (1976), 2213{2219. [78] M. J. Tindall, P. K. Maini, S. L. Porter, and J. P. Armitage. Overview of mathematical approaches used to model bacterial chemotaxis II: bacterial pop- ulations. Bulletin of Mathematical Biology 70 (2008), 1570{1607. [79] T. Vicsek, A. Czir ok, E. Ben-Jacob, I. Cohen, and O. Shochet. Novel type of phase transition in a system of self-driven particles. Physical Review Letters 6 (1995), 1226{1229. [80] N. Watari and R. G. Larson. The hydrodynamics of a run-and-tumble bac- terium propelled by polymorphic helical agella. Biophysical Journal 98(1) (2010), 12{17. [81] S. Webb, J. A. Sherratt, and R. G. Fish. Cells behaving badly: a theoreti- cal model for the Fas/FasL system in tumour immunology. Mathematical Bio- sciences 179 (2002) 113{129. [82] R. A. Weinberg. The Biology of Cancer. Garland Science, NY, 2007. [83] S. Wilson and D. Levy. A mathematical model of the enhancement of tumor vaccine e cacy by immunotherapy. Bulletin of Mathematical Biology (2012), accepted. [84] K. Yasui, M. Yamazaki, M. Miyabayashi, T. Tsuno, and A. Komiyama. Signal transduction pathway in human polymorphonuclear leukocytes for chemotaxis induced by a chemotactic factor. Distinct from the pathway for superoxide anion production. Journal of Immunology 152 (1994), 5922{5929. [85] R. Yu and D. Kaiser. Gliding motility and polarized slime secretion. Molecular Microbiology 63(2) (2007), 454{467. 143