APPLICATION OF C H E M O T A C T I C MODELS TO ALZHEIMER'S DISEASE by M A G D A L E N A L U C A B.Sc. (Mathematics) Transilvania University, Romania 1989 M.Sc. (Applied Mathematics) University of Manitoba, 1996 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF D O C T O R OF P H I L O S O P H Y in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Mathematics Institute of Applied Mathematics We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A January 2002 © Magdalena Luca, 2002 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of The University of British Columbia Vancouver, Canada Abstract This thesis presents the analysis of several chemotactic models developed for pattern formation in biological systems. The project was motivated by the complex interactions between glial cells, cytokines and proteins that cause neuronal death and eventual de-mentia in Alzheimer's disease. I investigate conditions that lead to periodic patterns and aggregation of glial cells and chemicals, as observed in the senile plaques of Alzheimer's disease. I have examined a hierarchy of one dimensional models, starting with a model in which there is a single chemoattractant inflammatory agent and the microglia cells. Thus, I constructed a minimal partial differential equations model similar to the Keller-Segel model for the chemotaxis of microglia and their involvement in secretion and uptake of various inflammatory factors. Stationary solutions of the system of equations are investigated using a method developed by Schaaf. The linear stability analysis of the system leads to the aggregation condition which can be biologically interpreted. As a result, I suggest possible methods that might affect plaque formation in Alzheimer's disease. Numerical simulations help corroborate theoretical results. A different form of the kinetics function which is part of the Keller-Segel model is sug-gested. Thus, the limiting kinetics function in the chemoattractant equation incorporates production or removal of the chemical factor by the cells, depending on the model pa-rameters. Also, it limits the sharp increase in the cell density which is well known to occur in this type of a chemotactic model. Finally, a model for chemotaxis of cells in response to a chemoattractant and a chemore-pellent is developed to include a greater level of detail. The system is reduced in complex-ity using quasi-steady state approximations, and explicit solutions for chemical diffusion are found using the Green's function method. The model predicts that local attraction and long range repulsion create periodic pattern formation. Biological interpretations of the linear instability condition lead to suggestions for therapy interventions in Alzheimer's disease. Using real biological parameters available in the literature, I compare and discuss the applicability of this particular model to the pattern formation observed in Alzheimer's disease. n Table of Contents Abstract ii Table of Contents iii List of Tables vi List of Figures vii Acknowledgements ix Dedication x Chapter 1. Introduction 1 1 Brief Description of the Biology 1 2 The Problem, the Motivation and the Approach 3 3 Review of Previous Work 4 4 Thesis Outline 10 Chapter 2. Alzheimer's Disease 12 1 The Cells in Alzheimer's Disease 12 1.1 Astrocytes 13 1.2 Morphological Properties of Microglia 15 1.3 Microglia, the Brain's Macrophage 16 1.4 The Harmful Microglia ; 18 2 The Chemicals in Alzheimer's Disease 19 2.1 Cytokines 19 2.2 Amyloid Associated Proteins 21 3 The Events in Alzheimer's Disease 21 3.1 Amyloid Plaque Formation 22 3.2 The Amyloid Cascade Hypothesis 25 4 Conclusion 26 Chapter 3. The Model and its Mathematical Background 28 1 Modelling Goals 28 2 The General Chemotaxis Model 29 3 Model Assumptions and Equations 30 4 The Cell Equation 31 4.1 The Stationary Solution for Cell Density 32 4.2 Conservation of Cells 33 5 The Diffusion Equation 34 ii i 6 The Reaction-Diffusion Equation 35 6.1 The Method of Green's Function 35 7 Numerical Simulations 40 7.1 The Finite Difference Method for Initial Boundary Value Problems 40 7.2 The Diffusion Equation 43 8 Conclusion 44 Chapter 4. The Keller-Segel Chemotactic Model 46 1 Model Equations and Mathematical Tools 46 2 Simplifications and Reductions 48 2.1 Space-Independent System 48 2.2 Time-Independent System 49 3 Linear Stability Analysis of the Space-Independent System 49 4 Linear Stability Analysis of the Full System 50 5 Biological Interpretations of the Aggregation Condition 53 6 Stationary Solutions 55 6.1 The Cell Equation 55 6.2 The Chemical Equation: The Reduced System 55 6.3 Fixed Points of the Reduced System 56 6.4 The Phase Plane cw 57 7 Numerical Simulations 59 7.1 Derivation of Parameters 60 7.2 Graphs 63 8 Conclusion 65 Chapter 5. Extensions of the Keller-Segel Model 68 1 Different Kinetics Functions for the Chemical Equation 68 1.1 Conditions Imposed on the Kinetics Function g(m, c) 70 1.2 Dynamics of the Chemotactic System 71 1.3 Numerical Simulations 73 2 Multiple Kinetics Functions for the Chemotactic Model 79 3 Conclusion 82 Chapter 6. The Chemoattraction-Repulsion Model 84 1 Model Equations and Mathematical Tools 85 2 Nondimensionalization 86 3 The Method of Green's Function 87 4 The Shadow System 90 5 Linear Stability Analysis 91 5.1 The Full System 91 5.2 The Shadow System 93 5.3 The Bifurcation Function 94 iv 5.4 The Parameter Space Aa 97 6 Biological Interpretations of the Aggregation Condition 100 7 Numerical Simulations 102 8 Conclusion 106 Chapter 7. Applicability of Models to Alzheimer's Disease 107 1 Parameter Derivation 107 2 Applicability of the Chemoattraction-Repulsion Model 114 2.1 Model Parameters 115 2.2 Comparison with Biology 120 2.3 Limitations of the Model 122 Chapter 8. Conclusion 124 1 Concluding Remarks 124 2 Future Work 126 Bibliography 129 Appendix A . Mathematical Background 137 1 Derivation of the Transient Solution : 137 2 The Steady State Solution 141 3 Derivation of the Finite Difference Equations 141 Appendix B. The Keller-Segel Chemotaxis Model 145 1 The Jacobian Matrix 145 2 The Lambert W Function 147 Appendix C. Extensions of the Keller-Segel Model 149 1 Different Kinetics Functions for the Chemical Equation 149 1.1 Reduced Systems 149 1.2 Linear Stability Analysis of the Space Free System 150 1.3 Linear Stability Analysis of the Full System 150 2 Multiple Kinetics Function for the Chemotactic Model 153 2.1 Reduced Systems 153 2.2 Linear Stability Analysis of the Space Independent System 153 2.3 Linear Stability Analysis of the Full System 154 Appendix D. The Chemoattraction-Repulsion Model 156 1 Derivation of the Nondimensional System 156 2 The Method of Green's Function 157 3 The Shadow System 160 4 Linear Stability Analysis of the Shadow System 161 v List of Tables 6.1 Dimensionless parameter values for the chemoattraction-repulsion model 102 7.1 Molecular weight of proteins and cytokines 108 7.2 Effective diffusivity coefficients of various molecules 108 7.3 Cell kinetics parameters for different cell types 110 7.4 Receptor kinetics parameters for cytokines 110 7.5 Binding rates (k+), unbinding rates (&_) and dissociation rates (KD) for cytokine receptors I l l 7.6 Concentration of cytokines produced by microglia or astrocytes I l l 7.7 Production rates of cytokines produced by microglia or astrocytes 113 7.8 Decay rates of cytokines 113 7.9 Additional parameters for cell movement 114 7.10Parameter values for the chemoattraction-repulsion model Equations 6.1 116 7.11Dimensionless parameter values for the chemoattraction-repulsion model 120 7.12Comparison between theoretical results and data from the biological literature. 121 vi List of Figures 2.1 Astrocytes in the human brain 13 2.2 The functional states of microglia 17 2.3 Plaques in an Alzheimer's brain (post-mortem) 23 2.4 When B and 7 secretases cleave the amyloid precursor protein molecule, they release the /3-amyloid protein 24 2.5 A schematic diagram of proposed plaque progression 25 3.1 Green's function for a set of parameters: x' = 1, b = 2, D = 8 40 3.2 Time behaviour of the chemical concentration c(x2i,t) 44 4.1 Dispersion relation for the Keller-Segel model, for different values of a 53 4.2 The phase plane for the Keller-Segel model 58 4.3 Plot of the linear function c — fi/x 62 4.4 Plot of the dispersion Equation 4.20 for the Keller-Segel model 62 4.5 Microglia distribution in the Keller-Segel model initiated with the stable mode n — 14 : 64 4.6 Stationary pattern formation in the Keller-Segel model initiated with the unstable mode n = 10 65 4.7 Spatial distributions of microglia and chemical in the Keller-Segel model 66 5.1 Sketches of the functions s(c) and h(c) in g(m,c) = s(c) — mh(c) 69 5.2 Dispersion relation for the extended System 5.1 74 5.3 Microglia distribution initiated with the stable mode n — 10 75 5.4 Microglia distribution initiated with the unstable mode n — 9 76 5.5 Microglia distribution initiated with the "stable" mode n = 10 77 5.6 Spatial distributions of microglia and chemical using random initial conditions. 78 5.7 Dispersion relation for study Case 1 81 5.8 Dispersion relation for study Case 2 82 6.1 Green's function G(x, x') -. 89 6.2 Plot of the dimensional decaying Green's function G(x, 0) 90 6.3 Graphs of the functions / 1 and f2 for different values of A and a 95 6.4 The bifurcation function a* for A = 2, a = 1.5, and the line 1/A2 for different parameter sets 96 6.5 The parameter space Aa 97 6.6 The bifurcation function a* in different regions of the parameter space Aa 98 6.7 The kernel K(x) in different regions of the parameter space Aa 99 6.8 The bifurcation function a* for A = 1.5, a = 2, and the line y = 1/A2 for A2 = 100 and A2 = 400 103 6.9 Spatial distributions of microglia and chemicals for the nondimensional chemotactic model, using random initial conditions 104 vii 6.10Spatial distributions of microglia and chemicals for the nondimensional chemotactic model, using random initial conditions 105 7.1 The bifurcation function a* for A = 1.5, a = 2, and the line 1/A2 for A2 = 100 and A2 = 400 118 7.2 Stationary spatial distributions of microglia and chemicals for the nondimensional chemoattraction-repulsion model, using random initial conditions 121 8.1 A schematic diagram of how neurons, microglia and astrocytes interact with proteins and cytokines 128 viii Acknowledgements I would like to express my deepest gratitude to my supervisor, Leah Edelstein-Keshet, for her advice and constant encouragement throughout the years I have spent at U B C . The work on this thesis was carried out under her expert direction and supervision. I cannot emphasize enough how much I appreciate her personal support in times that were very difficult for me. Thanks and good luck to Alexandra and Ricardo who were always there for me in the last two years! We exchanged ideas and worked closely in order to make progress on the Alzheimer's disease project. Many, many thanks go to Athan to whom I am so much indebted. He spent endless hours on helping me with the technical aspects of my thesis. Athan also gave invaluable suggestions regarding my research. I thank my advisory committee, Lionel Harrison, Robert Israel, Yue Xian L i and Brian Wetton. Their suggestions on many biological and mathematical aspects of my thesis proved to be invaluable to the development of my work. Many thanks to Chris Shaw and Ji l l McEachern for their insights into the biology of neurodegenerative diseases. Also, I must thank Alex Mogilner for the consultations and advice related to my thesis. I thank the Institute of Applied Mathematics for providing such resources as the lab and the new lounge. Thanks to the director of the Institute of Applied Mathematics, Bernie Shizgal, and to Roman Baranowski, the research/IT manager, for reestablishing the lab as a wonderful place of work! Warm thoughts go to all those people with whom I shared the lab and the lounge, for their friendship, interesting discussions and technical support: the alumni, Michele, Jennifer, Marty, Greg, Dave; the veterans, Adriana, Natalia, Stefan; and the new comers, Fatemah, Amy, Alyssa, Clive, Margaret, Ulrike, Marek, Stan. This research was funded by MITACS under the "Biomedical Models of Cellular and Physiological Systems and Disease" project, and partly funded by "In Silico Biosciences". During initial phases of research, I was funded by a University of British Columbia Grad-uate Fellowship (UGF) and Science Council of British Columbia G R E A T scholarship, as well as NSERC grant support (to Leah Edelstein-Keshet). And most of all, I have to thank my family, especially my husband Florin, for their endless support, understanding and continuous encouragements, and to Cassandra, who brought me so much joy and happiness! ix In loving memory of my father, Professor Constantin Iliescu (1930-2000), who taught me how to love science, whose work, life and love inspired me throughout my studies, and will continue to do so for the rest of my life. Chapter 1 Introduction In this thesis I am concerned with mathematical modelling of chemotaxis systems moti-vated by Alzheimer's disease. This introductory chapter gives an overview of the disease and describes the key biological aspects that will be translated into mathematical equa-tions. A review of previous work, relevant to my mathematical modelling approach, is also included. I conclude with a characterization of each chapter in the thesis. 1 Brief Description of the Biology Alzheimer's disease is a devastating neurodegenerative disease affecting almost one out of ten individuals who reach the age of 60, and more than half of those reaching the age of 85 (Cowley 2000; Jones 2000). The disease represents the eighth leading cause of death in America, lasting anywhere from 8 to 20 years (Nash 2000). It is characterized by a progressive decline of cognitive functions, fatal loss of mental functions, and eventual death. These reasons, stimulated by an increase of the elderly population in the world, motivated research in this area in the last two decades. As a result, progress on therapy interventions to slow and reverse the damage caused by Alzheimer's disease is very recent. (Most significant research has been done after 1985.) The disease was first investigated by Alois Alzheimer at the beginning of the 20th century (Banati and Beyreuther 1995). It is now well established that the brains of Alzheimer's 1 disease sufferers are affected by abnormal foci called senile plaques which are particu-larly concentrated in the hippocampus (Dickson 1997; Itagaki et al. 1989). The plaques are complicated lesions composed of extracellular deposits of the ^-amyloid protein, degenerating neurons and other non-neuronal cells called glia (Dickson 1997). Also no-table are abnormal structures inside neurons (neurofibrillary tangles) consisting of twisted filaments of the tau protein (Nilsson et al. 1998). Alzheimer himself observed and described both the plaques and the tangles in 1906 (Selkoe 1994). What causes senile plaques to form? One hypothesis suggests that the disease starts with local accumulation of soluble /3-amyloid protein into diffuse plaques (Banati and Beyreuther 1995). The concentration of soluble /3-amyloid is much higher in Alzheimer's disease than in normal individuals and it correlates with the severity of the disease (McLean et al. 1999). In time, the /3-amyloid builds up to form insoluble fibrils which stick together to create dense plaques (Banati and Beyreuther 1995). Inflammation caused by activation of glial cells by the /3-amyloid protein constitutes the next step in the plaque formation (Nilsson et al. 1998). Microglia and astro-cytes are the glial cells activated in this process. In vitro experiments indicate that reactive microglia start moving in response to concentration gradients of amyloid and many other chemicals (Itagaki et al. 1989). This biased movement is called chemotaxis and motivates much of the modelling approach described in the thesis. Both glial cells produce cytokines and additional chemical factors that intensify the inflammation process (Mrak et al. 1995; Mrak et al. 2000) and cause neuronal stress, further production and aggregation of the /5-amyloid protein and, finally, senile plaque formation (Nilsson et al. 1998). As a result, a cycle of amplified local inflammation can take place (Griffin et al. 1998), leading to local toxicity, stress and death of the neurons. I capture these phenomena by a simplified version, in which a single diffusible (toxic) 2 chemical is produced and removed. Cytokines represent a class of proteins involved in cell-to-cell communications that regulate the immune system (McGeer and McGeer 1999). Secretion of cytokines such as interleukin-1/3 (IL-1/5), interleukin-6 (IL-6) (Fiala et al. 1998), tumor necrosis factor-a (TNF-a) (Sawada et al. 1989; Lee et al. 1993) and S100/3 (Mrak et al. 1995) is well established, and many of these inflammatory cytokines are found in elevated levels in neurodegenerative diseases such as Alzheimer's disease and Parkinson's disease (Mogi et al. 1994; Fill i t et al. 1991; Mrak et al. 2000). The chemotactic attraction of microglia to amyloid (Davis et al. 1992), to various cytokines and proteins (Yao et al. 1990), and the general motility of these cells (Nolte et al. 1996) has been studied extensively. There is also evidence that T N F - a can induce a chemotactic repulsive response in cells (Chicoine et al. 1995). 2 The Problem, the Motivation and the Approach In this thesis, my modelling goal is to address and test the following hypothesis: can chemotaxis of microglia cells and their interactions with chemical factors lead to the formation of localized plaques, as seen in Alzheimer's disease? Thus, I will focus my attention on the early events in the formation of diffuse plaques in Alzheimer's disease. Although little is known about the early stages and the mechanisms of plaque formation in Alzheimer's disease, the research I did was commissioned by a start-up company, "In Silico Biosciences", interested in modeling Alzheimer's disease. The work presented here is a theoretical and mathematical treatment, and a simplified version of work done by Leah Edelstein-Keshet and Athan Spiros. Their Java simulations incorporate more details about the biology of Alzheimer's disease. 3 Biological problems can be presented mathematically in various levels of generality. If I were to describe minute detail, the physiological system would be represented by a three dimensional network that takes into consideration the special geometry of various cells, the multitude of chemical and cell types, and the events that they are implicated in. In practice, it is useful to consider simplified models first, learn the general features, identify the important parameters, and then add details when feasible, on the assumption that these important characteristics are retained in the complex system. The previous paragraph describes one of the main features of my modelling approach in this thesis: I begin my study with a simple mathematical model and continue with more complex ones. In addition, three important steps are followed when developing and analyzing the models: I start with an in-depth study of the biology of Alzheimer's disease to identify the most important biological variables; I then translate the biological details into mathematical equations; finally, I interpret biologically our mathematical results, and, when possible, compare them with specific aspects from the real biology. The modelling approach described above, i.e. the development of simple mathematical models incorporating two to three equations, has the advantage of allowing me to easily analyze and understand the dynamics of the system. In addition, I will be able to investi-gate and interpret the role of each parameter in the system. Unfortunately, such models incorporate only a few biological underlying mechanisms of the disease. In contrast, mod-els which include many biological details give rise to complex networks that are difficult to analyze, and difficult to understand from simulations, due to a large parameter space. 3 Review of Previous Work The purpose of this review is to highlight the background literature of the models analyzed in subsequent chapters. As I already mentioned, I am interested in modelling chemotaxis of cells in response to chemical gradients that lead to pattern formation. Therefore, I 4 focus my attention on reviewing previous studies of pattern formation in chemotactic models of biological systems. The most important reaction-diffusion theory of morphogenesis (the development of pat-tern and form) was developed by Turing (1952). He constructed a system of partial differential equations for chemicals which react and diffuse leading to the formation of steady state heterogeneous spatial patterns of chemical concentrations. Turing stated that if the chemicals have different diffusion rates, then spatially inhomogeneous pat-terns can develop under appropriate conditions. More specifically, diffusive instability is caused when a peak of activator concentration formed at some location starts to spread with a diffusion coefficient Da, initiating new peaks. Another chemical, the "inhibitor", diffuses more rapidly with a diffusion coefficient Di so the region surrounding the initial peak contains enough inhibition to prevent further peaks of activation. Patterns will form provided the range of inhibition is longer than the range of activation. This phenomenon is also known as lateral inhibition. The reaction-diffusion theory initiated by Turing was elaborated by many researchers interested in pattern formation. Thus, Prigogine and Levefer (1968) proposed a hy-pothetical system called the Brusselator and investigated spatially non-uniform chemical patterns. Later, Gierer and Meinhardt (1972) proposed several specific activator-inhibitor systems consisting of two diffusing substances. The activator enhances its own synthe-sis, as well as that of the inhibitor. The inhibitor makes both chemicals to decline in concentration. The interactions between these two substances can give rise to pattern formation under certain conditions. The Gierer-Meinhardt system was applied to the study of a variety of biological systems such as the modelling of Phyllotaxis, the arrange-ment of leaves on the stem of a plant; modelling of stomata, irregularly spaced pores on leaf surfaces; modelling of stripes in zebras and in the visual cortex (Meinhardt 1978). 5 More recently, Murray (1993) published an extensive study of pattern formation sys-tems. In his book "Mathematical Biology", he analyzes pattern formation determined by reaction-diffusion systems, by mechano-chemical interactions or by chemotaxis. I also mention here the book by Harrison (1993) which presents the kinetic theory of pattern formation in living organisms. The author discusses pattern formation from a philosoph-ical point of view, and also investigates the most important mathematical systems and their applications to pattern formation in biological systems. I continue this review by introducing the classic Keller and Segel (1970) chemotaxis model which describes the interactions between cellular slime molds and an attractant chemical that the cells secrete, cyclic adenosine monophosphate (cAMP). The slime molds are chemotactic to concentration gradients of c A M P and they aggregate, spontaneously forming "centers". The authors constructed a two dimensional model that takes into account the random motion and the chemotactic attraction of the cells (which assumed functional forms), the constant diffusion and the linear production and decay of cAMP. Keller and Segel (1970) discussed stability of the uniform steady state and suggested, for the first time, that the onset of cell aggregation can be regarded as an instability. I study this model extensively in Chapter 4. Lauffenburger and Kennedy (1983) analyzed a one dimensional chemotactic model for the inflammatory response of the polymorphonuclear (PMN) leukocytes (cells of the immune system) to bacterial invasion of tissue. The existence of both uniform and non-uniform steady states was proven. Using linear stability analysis, they showed that a leukocyte chemotactic response smaller than a critical value led to a non-uniform steady state, while a response greater than the same critical value stabilized the uniform steady state. This result of the study contrasted with the role of chemotaxis in aggregation of slime mold obtained by Keller and Segel (1970), since P M N leukocytes remove rather than produce the attractant. The authors also showed numerical simulations using parameter 6 values estimated from the experimental literature in order to evaluate whether their model predictions were biologically relevant. Tranquillo et al. (1988) also used a one dimensional model to measure cell migration pa-rameters, in particular, the random motility coefficient, //, and the chemotaxis coefficient, X- In the model, human P M N leukocytes were attracted to a small peptide, F N L L P . They proposed that x depends on attractant concentration. They compared their theoretical results with values obtained from experimental measurements and found good agreement. In another important paper, Rivero et al. (1989) used experimental observations of cell movement properties for two classes of cells: flagellar bacteria and P M N leukocytes. A generalized equation for the migration of chemotactic cell populations was developed from a probabilistic one dimensional model. The resulting population flux equations had different forms, and they were applied to two distinct experimental data sets for each of these cell types. Their main contribution was to relate cell population migration to de-tails of individual cell behaviour, since the Keller-Segel equations were not derived from an appropriate fundamental description of cell movement. Thus, population transport coefficients such as fi, the random motility coefficient, and x, the chemotaxis coefficient, were estimated from the models. The authors concluded that, under certain limiting con-ditions, the Keller-Segel equations represent reasonable approximations to their models. A simple cell-chemotaxis model for the generation of spatial patterns in cell aggregations was studied by Grindrod et al. (1989). Using no flux boundary conditions and varying the total number of cells in the model, the authors analyzed the local and global bifurcation of spatially heterogeneous patterns away from the homogeneous steady states. They also discussed the existence of periodic spatially distributed solutions for the cells and attractant in an infinite domain. A two dimensional chemotaxis model for spatial pattern formation was investigated ana-7 lytically and numerically by Maini et al. (1991), as a model for animal skin patterns. A logistic type growth for the cells and a Michaelis-Menten production rate of attractant were assumed. The authors studied steady states and bifurcating solutions, spatially heterogeneous solutions and the effects of domain size, aspect ratio and growth on the resulting patterns, showing that a simple chemotactic model could lead to the formation of a wide variety of patterns (stripes or spots). An excellent study of the dynamics of chemotactic systems was published by Othmer and Stevens (1997). The authors derived several general systems of partial differential equations that depend on how the movements of cells are affected by the local chemoat-tractant concentration. The systems were modelled by considering the cells to be random walkers producing a nondiffusible chemical which changes the environment of other cells. They investigated the dynamics of the systems and found the following: whenever the decay rate of the chemical was set to zero, collapse of the systems occurred; for nonzero values of the decay rate of the chemical, cells either aggregated or were concentrated in one sharp peak that determined blow up of the system. The models were applied to the study of myxobacteria movement in response to an attractant. Myerscough et al. (1998) proposed a cell-chemotaxis model to examine the role of the boundary conditions on spatial pattern formation. Most chemotaxis models use no flux boundary conditions, i.e. cells and attractant do not leave or enter the domain. Here the authors looked at the effects of mixed boundary conditions on the dynamics of pattern formation. Thus, numerical simulations for one dimensional and two dimensional mod-els were generated. The results showed that more general boundary conditions led to significant differences in the spatiotemporal properties of patterns compared with those exhibited under no flux boundary conditions. The differences were: no peaks forming on the boundary and the formation of stable patterns with asymmetric peaks. The authors applied their model to the spatiotemporal sequence of skeletal development in the chick 8 limb. Another study was written by Painter et al. (1999). A generalized Turing model with chemotaxis was developed to account for cell growth and movement. The authors an-alyzed the effects of these two processes on pigment patterning, and they applied the model to explain features of pattern formation in a growing system of angelfish called Pomacanthus. In addition, chemotaxis in response to chemical gradients led to aggre-gation of one type of pigment cell into a striped spatial pattern. Painter et al. (1999) argued that the applicability of classical' Turing models to biological pattern formation is limited because of the sensitivity of patterns to model parameter values. Most chemotactic models are developed to reflect the interactions between one cell type and one chemical. However, cells can respond chemotactically to several chemicals. This process is the focus of an investigation done by Painter et al. (2000). The model was a generalization of a previous one developed by Othmer and Stevens (1997). Numerical simulations for a one dimensional and a two dimensional model were generated. The authors showed that the presence of an attractant and a repellent in the system led to the formation of more complex spatial patterns. Thus, there were two types of patterns observed, depending on the parameter values: spots and thin or thick stripes. The authors applied and compared their model to the spots of a jaguar and the stripes of a lion-fish. A more recent paper written by Lee et al. (2001) presented the study of local and non-local mathematical models for biological phenomena. The authors constructed a one dimensional model for chemotaxis of myxobacteria, similar to Keller-Segel model. The mathematical tools they used for analysis were linear stability analysis, dimensional analysis, and perturbation theory. The main contribution of the study by Lee et al. (2001) lies in their investigation of non-local models. They showed that, when rate 9 differences between fast and slow local dynamics are large enough, local models are well approximated by non-local ones. They applied their models to two different biological systems, one from cellular biology and one from animal behaviour. 4 Thesis Outline The thesis contains six major chapters, in addition to the introductory chapter and the final one that includes concluding remarks and directions for future work. Chapter 2 gives a more detailed description of Alzheimer's disease than provided in this overview. Hypotheses that are relevant to my modelling purposes are the focus of my presentation. In chapter 3, I present the mathematical setting of the problem, i.e. the chemotactic system of partial differential equations. I show how the system is derived. I study subsets of the system, and indicate the mathematical tools used to analyze the models. Chapter 4 introduces the Keller-Segel chemotactic model, the simplest model developed for pattern formation in biological systems. I discuss stationary solutions, linear stability analysis, numerical simulations and possible therapy interventions in Alzheimer's disease. Chapter 5 explores different extensions of the previously introduced Keller-Segel model. In particular, I am interested in incorporating additional biological features into the system of equations. The more complex and novel chemoattraction-repulsion system of equations is analyzed in chapter 6. Here, I broaden my model to include more interactions between the cells and chemicals found in Alzheimer's disease. I study reduced systems derived from the original one, linear stability analysis and numerical simulations. I discuss the biological implications of the mathematical results. 10 Finally, chapter 7 connects parameters and features of the chemoattraction-repulsion model to the biology of Alzheimer's disease. Extensive research was conducted in or-der to identify real, biological parameter values and adapt them to model parameters. Some of the values were found directly in the medical literature, others were calculated using biological data, while the remaining unknown ones were estimated based on the mathematical results. I compare my model predictions with the real biology, and discuss limitations of the model. 11 Chapter 2 Alzheimer's Disease In the introductory chapter, I gave a brief description of Alzheimer's disease. In this chapter, I continue my analysis with a detailed presentation of aspects of the disease, focusing my attention on the features that will be mathematically modelled in subsequent chapters. Thus, I introduce the biological properties of the cells and chemicals that play an important role in Alzheimer's disease, and their interactions. I emphasize that my goal is not to make an exhaustive presentation of the pathology of Alzheimer's disease, but rather I would like the reader to have an understanding of some of the events leading to the disease. In addition, I will attempt to specify which aspects will be of direct relevance to my mathematical modelling. 1 The Cells in Alzheimer's Disease Nerve cell bodies and axons are surrounded by gl ia l cells. The name glia derives from the Greek word glue. Originally it was thought that the main purpose of glia in the central nervous system is holding the brain together (Kandel et al. 1991). As I will show, recent studies indicate that glial cells have other important functions, notably in the brain's immune system (Streit and Kincaid-Colton 1995b). There are between 10 to 50 times more glial cells than neurons in the central nervous system of vertebrates. More specifically, the human brain contains 100 to 130 billion 12 glia, 80,000 to 100,000 per cubic millimeter (Duffy 1983). Glial cells are divided into two major classes: macrogl ia and microglia . The macroglia subdivision consists of three predominant types of cells: astrocytes, oligodendrocytes and Schwann cells. I will be concerned mainly with astrocytes and microglia, as they play a central role in Alzheimer's disease. 1.1 Astrocytes Astrocytes are star-shaped cells that have the largest cell body and are the most nu-merous among glial cells (see Figure 2.1a). About 90% of all glial cells are astrocytes and oligodendrocytes (Duffy 1983). Astrocytes have several functions. The cells make intimate contact with neighbouring neurons and blood capillaries through their branches. These types of contacts form tight junctions also known as the b lood-bra in barrier that protects the brain by preventing many substances in the blood from entering the brain (Kandel et al. 1991). (a) (b) Figure 2.1: Astrocytes in the human brain, (a) Normal astrocyte. The cell body is indicated by an arrow. Scale bar 10 pm. [From Zhou et al. (1990), by copyright permission.] (b) Reactive astrocyte. The cell body is 15 to 25 fim. [From Duffy (1983), by copyright permission.] Astrocytes recycle certain neurotransmitter substances around neurons, notably the ex-citatory transmitter, glutamate. In addition, astrocytes absorb and buffer the excess 13 K+ ions released during periods of intense neuronal activity, thereby protecting nerve cells from receiving too much stimulation (Kandel et al. 1991). Astrocytes also release neuronal growth factors that promote neuron growth and survival (Giulian 1995). For many years it has been known that Alzheimer's disease is associated with astrocytic gliosis (Dickson et al. 1995). Gliosis is the term used to describe the response of reactive glial cells to neuronal damage (Streit et al. 1999). When lesions take place in the brain, astrocytes proliferate, become reactive (see Figure 2.1b), and migrate to the injury site (Zhou et al. 1990; Faber-Elman et al. 1995; Faber-Elman et al. 1996), where, along with microglia, remove neuronal debris and fill the spaces that were occupied by neurons and their processes (FitzGerald 1992). Experiments in vitro indicate that astrocyte migration is blocked by two cytokines, T N F - a and IL-1 (Faber-Elman et al. 1995; Faber-Elman et al. 1996). Astrocytes play important roles in Alzheimer's disease. Thus, the cells can be found in senile plaques that have a dense amyloid core, and are mostly located at the periphery of the plaque (Itagaki et al. 1989). In addition, astrocytes and microglia are functionally inter-related: one cell type affects the other. Most central nervous system inflamma-tory and degenerative disorders are associated with reactive microglia and astrocytes. As already mentioned, activated astrocytes release growth factors that further promote microglial proliferation and activation (Dickson 1997). Astrocytes have been shown to inhibit phagocytic properties of microglial cells, thus amplifying local toxicity (Smits et al. 2000). More importantly, astrocytes produce inflammatory cytokines such as IL-6 (Lee et al. 1993) and S100/3 (Mrak et al. 1995), which contribute to a local cycle of neurotoxicity. They also synthesize and release other harmful proteins (such as apolipoprotein E and ai-antichymotrypsin (Mrak et al. 2000)). 14 1.2 Morphological Properties of Microglia Microglia are ubiquitous throughout the brain of humans in both gray and white matter. They number about 10% of the total glial cell population, and therefore they are at least equal to the number of neurons, that is 1 0 u (Streit 1995a). Microglia cells are small in comparison with astrocytes. They have elongated or triangular nuclei, a cell body that measures approximately 10-15 pm (Streit 1995a) and extensive branching of cytoplasmic processes. These processes expand into every region of the brain and touch neurons and astrocytes, but not one another (Streit and Kincaid-Colton 1995b). Unlike other brain cells, microglia originate outside the central nervous system, i.e. in the bone marrow, and enter the circulation as monocytes1 during late fetal life (FitzGerald 1992). At any subsequent time, more microglia can be generated by blood monocytes or they may locally divide (McGeer and McGeer 1995). This proliferation usually occurs during a period of high phagocytic demand. Phagocytosis represents the process of en-gulfment and digestion of foreign particles by specialized cells, including macrophages. We will see that microglia, under certain conditions, become phagocytic cells within the brain. Microglia in the adult brain can take various morphological appearances that correlate with distinct functional states (Streit et al. 1999). As seen in Figure 2.2, there are at least three identifiable states: resting or ramified microglia found in the normal brain; activated or reactive microglia that can be found after neuronal injury occurs; and phagocytic microglia that function as the macrophages of the brain. In addition, as we age, microglia cells become increasingly reactive (Sheng et al. 1998), and they also begin to stick to each other to give rise to small microglial clusters (Streit et al. 1999). If such clustering of microglia increasingly progresses, it can culminate in the formation LMonocyte = A variety of white blood cell. 15 of senile plaques (Streit et al. 1999). I should emphasize that, throughout the thesis, I will be interested in the reactive functional state of microglia cells. For my modelling purposes, it is important to note that Alzheimer's disease patients have a greater density of microglia cells than normal individuals (Mackenzie et al. 1995). Thus, experiments conducted in vitro using post-mortem tissue from the temporal neocortex brain region, show that Alzheimer's disease patients had 0.5 cells per 104 \im2 compared to 0.3 cells per 104 /im2 for normal individuals displaying both diffuse and dense plaques (Mackenzie et al. 1995). Furthermore, the hippocampus of Alzheimer's disease patients shows a much higher density of reactive microglia cells (Itagaki et al. 1989): in vitro experiments indicate that there are 100 to 350 reactive microglia cells per 104 \irrv2. 1.3 Microglia, the Brain's Macrophage The major role microglia have during embryo development is to remove dead cells since the growing fetus generates more neurons and glial cells than it needs (Streit and Kincaid-Colton 1995b). In the normal adult brain, the most important function of microglial cells is to serve as the brain's macrophage cells by reacting to injury and disease. Resting mi-croglia are the first cells, followed by astrocytes, to respond almost instantly (within min-utes) to disturbances in their environment by moving toward the injured sites (Gehrmann et al. 1995). Chemotaxis is a complex process in which the direction of a moving cell is biased by the concentration gradient of some chemical substance (Wilkinson and Haston 1988). The chemical is a chemoattractant (or a chemorepellent) for the cell if the cell moves in the direction of (or away from) a higher concentration (Gradon and Podgorski 1995). Chemoattractants lead to positive chemotaxis, while chemorepellents lead to negative chemotaxis. The ability of a cell to change its direction of movement in response to a concentration gradient can be represented by the chemotactic coefficient x- The value 16 - 5 ^ i > 3 PEniNEUTOMAL MICROGLIA PHAGOCYTIC MICPOGUA Figure 2.2: The functional states of microglia. Resting microglia are highly rami-fied (top). Neuronal injury activates microglia. An immediate sign is enlargement of processes and change of shape. Depending on the struc-ture of their environment, microglia may take different appearances (middle). If neuronal death occurs, microglia become phagocytic and try to remove the debris (bottom). [From Streit and Kincaid-Colton (1995b), by copyright permission.] 17 of this coefficient depends on many cell properties: random motility, velocity, receptor kinetics between the cell and the chemical (Gradon and Podgorski 1995). I note that all my mathematical models take into consideration chemotaxis. In the brain, chemotaxis results in attraction of microglia to injured areas (Gilad and Gilad 1995; Streit and Kincaid-Cblton 1995b). The ramified processes of microglia enable them to be very sensitive to environment changes, and therefore they easily detect any neuronal injury. The cells then migrate to the injury site, assess the severity of damage, and possibly become phagocytic (Streit and Kincaid-Colton 1995b). Reactive microglia can revert to their resting state if the detected injury is mild or reversible. If the injury is severe and neuronal death is detected, microglia continue to retract their processes and become phagocytic, trying to clean up the site by removing dead neurons (Streit and Kincaid-Colton 1995b). I should note that microglia can become phagocytic in response to many harmful chem-icals. Ard et al. (1996) showed, for example, that microglia are capable of rapidly and effectively removing soluble ^-amyloid protein from tissue and /3-amyloid aggregates formed in diffuse plaques. Furthermore, Shaffer et al. (1995) indicated that by clearing soluble /3-amyloid, fibrils of the protein cannot form, and therefore microglia play an im-portant role in the defense against the formation of dense plaques containing /3-amyloid. The removal of /3-amyloid by microglia cells will be later considered in my mathematical models. 1.4 The Harmful Microglia In contrast, highly activated microglia can have serious bystander effects by which healthy cells are harmed: microglia release many substances that, at high levels, are dangerous or toxic to neurons (Streit and Kincaid-Colton 1995b). Microglia cells, in response to the stimulating ^-amyloid protein, secrete cytokines, free radicals of oxygen 18 and nitric oxide, excitatory amino acids like glutamate, which all damage neurons and lead to a cascade of inflammatory events (Barron 1995; Banati et al. 1993; Smits et al. 2000). Several studies (Mackenzie et al. 1995; Stalder et al. 1999) show that, although neu-rons are the main source of the /3-amyloid protein, microglia are also able to produce /3-amyloid. The capability of microglia to secrete /3-amyloid protein represents a contro-versial hypothesis (Shaffer et al. 1995). 2 The Chemicals in Alzheimer's Disease The local inflammation that occurs in the Alzheimer's brain is due to the ability of endogenous cells, such as neurons, astrocytes and microglia, to synthesize many pro-inflammatory chemicals involved in Alzheimer's disease lesions (McGeer and McGeer 1999). For my modelling purposes, I am interested in two types of chemicals produced by these cells: cytokines and amyloid associated proteins. 2.1 Cytokines Cytokines play an important role in the initiation, propagation, regulation and suppres-sion of immune and inflammatory responses in the central nervous system (Benveniste 1995). They are low molecular weight proteins that can be produced by, and have mul-tiple biological effects on different cell types (Benveniste^ 1995). In the normal brain, cytokines are present at very low concentrations and are rapidly induced when injury occurs. The level of concentration in tissue is indicative of the type and intensity of inflammatory activity. Cytokines are removed from the brain due to uptake by cells through receptors2 found on the cell's surface. 2Receptor = Surface-bound protein molecule that function to acquire information from the extracellular space and transmit this information into the cell through the plasma membrane. 19 The major classes of cytokines are the interferons, interleukins, tumor necrosis factors and colony stimulating factors (Mrak et al. 2000). For my modelling purposes, I am mainly concerned with interleukin-1, interleukin-6, and tumor necrosis factor-a. I mention here that the molecular weight of proteins is given in kDa which represents the atomic mass unit of a protein and it is equal to 1.67 x 1 0 - 2 4 kg. Interleukin-1 (IL-1) has a molecular weight of 17 kDa and it is a cytokine produced commonly by activated macrophages in the body and by microglia in the brain (Ben-veniste 1995). IL-1 is expressed in two major forms, I L - l a and IL-1/3 (Benveniste 1995). IL-1 mediates a variety of processes in response to inflammatory disease. Thus, it plays an important role in the brain since it can induce in many cells of the central nervous system the production of various other cytokines, such as IL-6, T N F - a and IL-1 itself (Nilsson et al. 1998). IL-1/3 is known to stimulate astroglial proliferation in response to brain injury (Benveniste 1995). More importantly, IL-1 induces astrocytes and neurons to produce more /3-amyloid which leads to the formation of amyloid fibrils in Alzheimer's disease (Nilsson et al. 1998). Interleukin-6 (IL-6) is a multifunctional cytokine involved in the regulation of inflam-matory responses in Alzheimer's disease (Nilsson et al. 1998). IL-6, a 26 kDa molecule, is secreted by a variety of activated cells, including microglia and astrocytes (Benveniste 1995). In the brain, astrocytes are the major producers of IL-6 when they are stimulated by IL-1/3 (Lee et al. 1993). Elevated levels of IL-6 are observed when senile plaques start to form in the brain of patients with Alzheimer's disease (Nilsson et al. 1998). IL-6 increases blood-brain barrier permeability, induces gliosis, and stimulates production of pro-inflammatory cytokines (Merrill and Benveniste 1996). Tumor necrosis factor-a (TNF-a) is a major pro-inflammatory cytokine found in the central nervous system. This 17 kDa molecule can be secreted by stimulated astrocytes 20 and microglia (Benveniste 1995). T N F - a can affect immune reactions and it can stimulate many cell types to produce cytokines, including IL-1, IL-6 and T N F - a itself (Benveniste 1995). T N F - a also has the ability to induce signals that trigger programmed cell death (called apoptosis) (Venters et al. 2000). In addition, when the brain is exposed to IL-1 or T N F - a , the blood-brain barrier permeability is altered and more inflammatory cells migrate into the central nervous system, leading to local neurotoxicity (Lee et al. 1993). 2.2 Amyloid Associated Proteins Amyloid deposits contain, in addition to /3-amyloid peptides, other proteins that are nor-mally observed during inflammation. There are at least two types of amyloid associated proteins: apolipoprotein E and ax-antichymotrypsin (Nilsson et al. 1998). A p o l i p o p r o t e i n E (ApoE) is a 37 kDa lipoprotein 3 primarily produced by astrocytes and found in senile plaques (Krul and Tang 1992). The ApoE4 genotype is associated with increased risk for late onset familial and sporadic Alzheimer's disease (Dickson 1997) . In addition, ApoE4 is active in promoting the aggregation of diffuse amyloid into fibrillar amyloid (Mrak et al. 2000). It also contributes to the insolubility of extracellular /3-amyloid protein deposits and thus inhibits its clearance (Yankner 1996). a i -an t ichymotryps in (ACT) is an integral element of amyloid filaments (Nilsson et al. 1998) . The protein is produced by astrocytes surrounding amyloid plaques and its pro-duction is induced by IL-1 (McGeer and McGeer 1995). 3 The Events in Alzheimer's Disease The progression of neurodegenerative changes that take place in the brains of adults with Alzheimer's disease are not well understood. The diagnosis of the disease is difficult to l ipoprotein = A class of proteins found in blood plasma that are combined with other fats or lipids. 21 make and it is done through physical and psychological exams, but it can be confirmed only following examination of brain tissue after death. What makes the study of the events in Alzheimer's disease even more difficult is that there are no Alzheimer's disease conditions in experimental animals (McGeer and McGeer 1999). Mice do show some similar features, but it is hard to make the connection to humans (McGeer and McGeer 1999). In the series of events leading to Alzheimer's disease, Griffin et al. (1998) suggests that IL-1 represents one of the initial key players in a cascade of incidents. IL-1 facilitates neuronal synthesis and processing of the amyloid precursor protein that promotes depo-sition of /3-amyloid; it activates glial cells, in particular microglia and astrocytes, which in turn contribute to release of more inflammatory cytokines and amyloid associated pro-teins. Thus, overexpression of IL-1 produces feedback amplification of this inflammatory process that culminate with Alzheimer's disease. In this section, I attempt to explain what are some of the important events that lead to Alzheimer's disease. 3.1 Amyloid Plaque Formation A m y l o i d plaques are the major biological markers of Alzheimer's disease (Figure 2.3b), and their formation represent the focus of my mathematical investigation. Figure 2.3a shows the distribution of such amyloid plaques in the hippocampus of an Alzheimer's disease patient (the biological details observed in this figure will be later compared with my mathematical predictions). The plaques range in size from approximately 10 to 100 fim, while the distance between them is about 50 to 200 \xm. There are approximately 0 to 30 amyloid plaques per 104 pm2 of tissue (Itagaki et al. 1989). The amyloid plaques found in the post-mortem brains of Alzheimer's disease patients consist of /3-amyloid protein (Goedert 1993). Alzheimer described the protein as a "peculiar substance", but it was named and characterized only in the 80s (Selkoe 1991). 22 (a) (b) Figure 2.3: Plaques in an Alzheimer's brain (post-mortem), (a) Diffuse (lighter grey) and dense (darker grey) plaques. Bar 55 fim. (b) Single amy-loid plaque surrounded by astrocytes. [From Itagaki et al. (1989), by copyright permission.] The /3-amyloid protein is an approximately 38-42 long amino acid protein with a molec-ular weight of 4.5 kDa (Dickson et al. 1995). The protein is a collection of peptides derived from the consecutive cleavage of the amylo id precursor pro te in at two sites. The amyloid precursor protein is a normal protein produced by healthy neurons (Cowley 2000). There are three proteases4, known as the a, 8 and 7 secretases, that can split the amyloid precursor protein (Mesulam 1999). When the 8 and 7 secretases cleave the amyloid precursor protein, they release the /3-amyloid protein (see Figure 2.4). The /3-amyloid protein in Alzheimer's disease seems to be both of neuronal and non-neuronal origin (Wisniewski and Wegiel 1995). Studies suggest that the neuron is the source of amyloid in diffuse plaques, while microglia cells produce amyloid in dense plaques (Mackenzie et al. 1995; Stalder et al. 1999). A l l individuals produce the /3-amyloid protein, but it is not known what purpose it serves. Unfortunately, under certain circumstances, it can build up in the extracellular spaces surrounding neurons Protease = A digestive enzyme that causes the breakdown of protein. 23 Figure 2.4: When 8 and 7 secretases cleave the amyloid precursor protein molecule, they release the /3-amyloid protein. to form plaques (Cowley 2000). A typical Alzheimer's disease patient may produce /3-amyloid at the same rate as a healthy person. The difference between the two individ-uals stems from how each is able to dispose of the /3-amyloid. The evolution of an amyloid plaque is shown in Figure 2.5, based on a schematic rep-resentation of plaque progression proposed by Griffin et al. (1997). Initially, we have diffuse deposits of soluble /3-amyloid that contain no dead neurons (difuse non-neuritc plaques). Later, the /3-amyloid protein is condensed and abnormal neurons appear to be part of such a deposit, giving rise to difuse neuritic plaques. Inside these the /3-amyloid protein slowly transforms into insoluble fibrils and forms a compact core of /3-amyloid, leading to the formation of a dense core neuritic plaque. With further progression, there is a loss of neuritic elements and formation of a dense core of /3-amyloid called a dense core non-neuritc plaque. These dense core non-neuritic plaques are considered to represent the final state of plaque evolution. As already mentioned, amyloid plaques also contain reactive glial cells (Griffin et al. 1998). Reactive astrocytes are usually located at the margins of amyloid deposits, while reactive microglia are situated in the center of the amyloid deposit (Itagaki et al. 1989). 24 Ditiuse Dilluse Dense core Dense core n o n - i K u r i t i c neuiilic ncuriiic non-neurilic plaque plaque plaque plaque Figure 2.5: A schematic diagram of proposed plaque progression. Cytokines and amyloid associated proteins released by reactive microglia and astrocytes accelerate the /3-amyloid aggregation by transforming and stabilizing the amyloid fibrils to form cores of dense plaques (Nilsson et al. 1998). However, amyloid plaque formation depends not only on the rate of /3-amyloid production, but also on the rate of its removal (Shaffer et al. 1995). Overproduction of the peptide could overcome the ability of microglia to remove soluble or fibrillar /3-amyloid. Furthermore, astrocytes might inhibit removal of /3-amyloid protein by microglia (Shaffer et al. 1995; Smits et al. 2000). The hypothesis that /3-amyloid deposition is the central event in Alzheimer's disease has been challenged on the basis of studies that show that the number of neurofibril-lary tangles correlate with the severity of dementia more accurately than plaque density (Yankner 1996; Hardy 1997). It is also possible that /3-amyloid may not be sufficient to cause Alzheimer's disease, but rather may act synergistically with other age-related pathological factors such as oxidative stress (Lockhart et al. 1994). 3.2 The Amyloid Cascade Hypothesis The leading model that explains the causes of developing Alzheimer's disease is the amylo id cascade hypothesis (Nilsson et al. 1998; Smits et al. 2000). The hypothesis asserts that overproduction or insufficient clearance of /3-amyloid sustains the deposition and formation of amyloid plaques as the central event in Alzheimer's disease. Further, 25 according to this hypothesis, the /3-amyloid is responsible for producing all other features of the disease; it causes neuronal degeneration and, ultimately, dementia. Nilsson et al. (1998) describe the three steps of the amyloid cascade hypothesis: accumu-lation of /3-amyloid into diffuse plaques; inflammation caused by microglial and astrocytic activation; production and release of inflammatory cytokines (IL-1, IL-6, TNF-a) and amyloid associated proteins (ApoE, A C T ) . These inflammatory cytokines and proteins' intensify neuronal stress that enhances /3-amyloid protein production and aggregation, and, finally, senile plaque formation (Nilsson et al. 1998; Smits et al. 2000). Thus, a feedback amplification and propagation of this cascade of events is produced, leading to local neurotoxicity. I should mention that there are several other competing hypotheses for the causes of Alzheimer's disease. For example, it is known that neurobrillary tangles distributions in the brain correlate better with the severity of Alzheimer's disease (Mesulam 1999), and therefore much attention should be given to their study. The age dependence of Alzheimer's disease has given rise to speculation that oxidative stress may play a role in the pathology of the disease (Yankner 1996). The generation of oxygen free radicals has been implicated in other neurodegenerative disorders such as Parkinosn's disease and amyotrophic lateral sclerosis (Yankner 1996). Another important hypothesis proposes that aged individuals with an AP0E4 genotype may have accelerated loss of neuronal functions and a significant increase in the number of amyloid plaques. In addition, FitzGerald (1992) indicates that severe head injuries cause inflammation in the brain which can lead to Alzheimer's disease. 4 Conclusion In this chapter, I described in detail one of the hypotheses that leads to Alzheimer's disease, in particular the amyloid cascade hypothesis . Clearly, the progression of events 26 is not yet fully understood. But the disease has specific features such as deposition of /3-amyloid protein, plaque formation, reactive glial cells and elevated concentrations of cytokines and amyloid associated proteins. In the next chapters, I will attempt to mathematically model the initial events of diffuse plaque formation: chemotaxis of mi-croglia in response to /3-amyloid and cytokines; removal and production of /3-amyloid and cytokines; inflammation cycle; microglia aggregation. 27 Chapter 3 The Model and its Mathematical Background In the previous chapter, I described the biological aspects of my topic of interest, i.e. mathematical modelling of Alzheimer's disease. In this chapter, I introduce the math-ematical setting for my problem. I present a detailed analysis of various subsets of the complex system of partial differential equations. Several mathematical tools that will be used throughout the thesis are also introduced here. 1 Modelling Goals The progression of Alzheimer's disease is complex and, as yet, poorly understood, and therefore I will not aim to address this full complexity. I will be primarily interested in the early events in the initiation of diffuse plaques. My goal is to understand the role of microglia (represented by its density m), their chemotaxis, uptake and secretion of chemical factors. I am particularly interested in the cycle of inflammation, in which local toxic chemical results in further production of toxic chemical and in more local microglial clustering. With this in mind, I omit the later stages of plaque development in which astrocytes are involved. I will also avoid a full representation of all known neuroprotective and neurotoxic influences, and represent these in a single variable for local toxicity, c. This variable could represent the concentration(s) of amyloid, cytokines, or some combination. Elevated levels of c are here associated with local neuronal death. It is important to note 28 that the chemical c is a chemoattractant for the microglia cells. Thus, I ask: what types of interactions between microglia and the toxic chemical factors would lead to the formation of localized aggregates? 2 The General Chemotaxis Model The equations that I use in this thesis to mathematically model the cell interactions are classic, being the traditional Keller-Segel equations for chemotaxis (Keller and Segel 1970,1971a,1971b). Almost 30 years after its proposal, Keller and Segel's model remains by far the most popular model for pattern formation as a result of chemotaxis. The model involves the cell density, m(x,t), and the chemoattractant c(x,t), where x € R 3 and t are the spatial coordinate and time respectively, and consists of two coupled reaction diffusion equations with the following general form: ^ - - V - J m , (3.1a) dc — = - V - 3c + g(m,c) , (3.1b) where J m is the flux of cells, J c is the flux of chemicals, and g(m, c) represents the net chemical production, dm/dt and dc/dt denote the change of m and c per unit time, while V - represents the divergence of a vector and leads to a scalar quantity. Detailed deriva-tions of this model are given by Edelstein-Keshet (1988) and Keller and Segel (1971a). In Equation 3.1a, the flux of cells J m stems from two different mechanisms: first, the cells move randomly throughout the domain. Second, the cells undergo chemotaxis, i.e direc-tional movement of cells up a gradient in chemical concentration. Thus, it is assumed that there are two contributions to the flux term J m , namely random Fickian motility with J m o t i i i t y = ~^{m^ c)Vm, where p(m, c) represents the random motility coefficient of microglia, and chemotactic type flux with J chemotaxis = ~x{mi c)mVc, where x(m, c) is the chemotactic coefficient of microglia. For the chemical Equation 3.1b, I assume that the only flux is Fickian diffusion: Jdiffusion = —D(m, c) Vc , where D(m, c) is the diffusion 29 coefficient of the chemical Using these assumptions, System 3.1 becomes random cell movement chemotaxis dm ~dt dc dt V - ( /i(m, c ) V m - x(m, c ) m V c ) V • (D(m,c)Vc) + g(m,c) . (3.2b) (3.2a) chemical diffusion kinetics function As already mentioned, I am interested in plaque formation in the brains of Alzheimer's disease patients. The plaques usually form in specific parts of the brain, for example in the hippocampus. Therefore, I consider Equations 3.2 on a finite domain T with Neumann (no flux) boundary conditions, namely where n is the unit outward normal to the boundary dT. With the specification of the boundary conditions, I conclude my theoretical introduction of the general chemotaxis model. Thus, the model is governed by a system of equations, as given by 3.2, coupled with the boundary conditions 3.3. In the next section, I will adapt this to a one dimensional chemotaxis model to be studied in subsequent chapters. 3 Model Assumptions and Equations For simplicity, I consider a one dimensional setting, 0 < x < L, where L is a typical dimension in the brain. I define rn(x, t) to be the density of activated microglia at time t and c(x, t) to be the concentration of the toxic chemical along my domain. The chemical c(x, t) plays the role of chemoatrractant for microglia. I assume that microglia cells undergo random motion with constant random motility coefficient fi and superimposed on chemotaxis towards high chemical concentration with constant chemotactic coefficient x- No proliferation of cells is considered. The chemical factor can diffuse in the region with constant diffusion coefficient D. The term of interest n • V m = n • V c = 0 , for x € dT , (3-3) 30 is the rate of uptake/secretion of the chemical factor by the cells, represented by a kinetics function g(m,c). Thus, in a one dimensional setting, with constant coefficients, System 3.2 becomes dm d2m d ( dc\ . . tW = ^ - x t e \ m d i ) ' (3-4a) dc d2c at = + »("•.=)• ("") The one dimensional no flux boundary conditions follow from the general boundary conditions 3.3: dm dx dc dx = 0 , (3.5a) x=0,L = 0 , (3.5b) x=0,L so that cells and chemical do not leave or enter the domain. I will show that this implies that the total number of cells in the domain is conserved. In addition, I provide the following initial conditions: m(x, 0) = u(x) , 0 < x < L , (3.6a) c(x, 0) = w{x) , 0 < x < L , (3.6b) where u(x) and w(x) are some sufficiently smooth functions. 4 The Cell Equation In order to analyze the behaviour of the partial differential equations that form Sys-tem 3.4, it is important to understand the dynamics of subsets of this system. In this section, I analyze the cell Equation 3.4a. This equation has no known analytical solution. However, it is possible, and insightful, to determine its stationary solutions, as shown below. 31 4 .1 The Stationary Solution for Cell Density I define the stationary solutions, or steady state solutions, of a system as the values of the variables for which there are no temporal changes in the system (Edelstein-Keshet 1988). In general, for any given equation, they are obtained by setting the time derivative to zero. Therefore, my stationary cell solutions can be calculated by letting dm/dt = 0 in Equation 3.4a. Suppose that the chemical concentration is given by some predetermined function of the spatial variable c = c(x) which is time independent. I have Integrating once leads to: where K\ is an arbitrary constant of integration. The expression above represents flux of cells, and using the Neumann boundary conditions 3.5,1 find that K\ = 0 so that the equation of interest is ^ - x m Y x = ^ ( 3 ' 9 ) Assuming that m ^ 0 for x G [0, L], I can separate variables in Equation 3.9 to obtain 1 dm x 9c m dx ii dx Integrating both sides leads to (3.10) ln(m) = - c + K2 , (3.11) A4 where K2 is an arbitrary constant of integration. After renaming the integration con-stants, the above can be rewritten in the form m = Meic, (3.12) 32 where M > 0 is constant and both m and c are functions of the spatial variable, i.e. m(x) = M e i c { x ) . (3.13) The cell Equation 3.13 is the non-uniform stationary solution, or non-homogeneous stationary solution, linking cell density to chemical concentration. Biologically, such solutions represent a non-uniform distribution of cells, reached in the limit after a long time, in my domain. In contrast, spatially uniform stationary solutions, or ho-mogeneous stationary solutions, are constant over space and time. These solutions characterize a uniform distribution of cells throughout the domain which does not change in time. Note that the cell density depends exponentially on the chemical concentration. This exponential relationship that exists between the stationary cell distribution and station-ary chemical concentration, as given by Equation 3.13, is an important feature of the chemotactic model and it will be addressed often in the following chapters. 4.2 Conservation of Cells As mentioned before, I show that Equation 3.4a, together with the boundary conditions 3.5, implies conservation of cells, i.e. of the total amount of m in the domain. Integrating equation 3.4a over the domain 0 < x < L , I obtain: dt . / n dx dm dc ox dx dx . 'o u o JO. Interchanging the order of integration and differentiation on the left-hand side, and using the Fundamental Theorem of Calculus on the right hand side leads to: L r o c - -\ \L 77- / mdx dtj0 dm dc dx dx Defining the total number of microglia in the domain, M(t) = m(x,t) dx, and using the Neumann boundary conditions 3.5, I find that d M —r- = 0 dt 33 This verifies that the number of microglia in the domain are conserved by the system. 5 The Diffusion Equation I consider simple chemical diffusion over a one dimensional domain with no flux boundary conditions, i.e. I look at the following equations: dc d2c di dx^ ' dc dx 0, (3.14) (3.15) x=0,L with initial condition c(x, 0) = w(x) , 0 < x < L . (3.16) The problem described by Equations 3.14, 3.15 and 3.16 is called in the literature (Za-uderer 1989; Boyce and DiPrima 1992) an in i t i a l boundary value problem. For this specific initial boundary value problem, I look for solutions of the form c(x, t) = c(x) + CT(X, t) , (3.17) where c(x) is the steady state solution of Equations 3.14 and 3.15, and cr{x,t) is called the transient function. Using the method of separation of variables as described in Boyce and DiPrima (1992), I find, as outlined in Appendix A l , that the solution to my initial boundary value problem is given by oo c(x, t) = c + ^ cn cos(\nx)e~DXnt , (3.18) n=l with 2 fL = — I w(x) cos(\nx) dx , n = l , 2 , 3 . L Jo (3.19) Here c = c(x) is the homogeneous steady state of the Equation 3.14, as found in Ap-pendix A2, and the second term in Equation 3.18 represents the transient function 34 cr(x,t). Xn = nir/L are called the eigenvalues of the initial boundary value problem, Xn(x) = cos(A nx) are the corresponding eigenfunctions, and n is called the mode. Thus, the study of the diffusion Equation 3.14 gives us information about how the chem-ical behaves when there is no production or removal of substances. In Equation 3.18, as t —> oo, CT(X, t) —> 0, and thus c(x, t) —> c. Therefore, in the absence of the kinetics function g(m, c), as time passes, the chemical always diffuses and reaches a homogeneous distribution c. In section 7.2, using numerical simulations, I will show this phenomenon graphically. In the next section, I study the diffusion of the chemicals in the presence of the kinetics function g(m, c). 6 The Reaction-Diffusion Equation In the previous section, I found the analytical solution to the diffusion Equation 3.14. In this section, I examine analysis of the reaction-diffusion equation | = D g + , ( m , c ) , (3.20) where g(m, c) has the specific form of a linear function g(m, c) = am — be. I will show that, given a specified time-independent distribution of cells m(x) (assumed static), the chemical distribution is prescribed by /•OO c(x) = / e-V&x-x'\m(x') dx' . (3.21) This form is obtained using the method of Green's function (Zwillinger 1989), as shown in the next section. A reader familiar with this technique may skip to Section 7. 6.1 The Method of Green's Function One of the mathematical tools that will be used later in the thesis is the method of Green's function. In order to solve Equation 3.20 using this method, a couple of 35 conditions must be satisfied. First, the chemical should have a stationary distribution, i.e. the left-hand side of Equation 3.20 is equal to zero: d2c 0 = £ > — -bc + am . (3.22) ox1 Another mandatory condition is that the kinetics function g(m;c) must be linear in c, condition which is already satisfied by the form of g(m, c). Equation 3.22 is a second or-der, linear, non-homogeneous ordinary differential equation in c. The method of Green's function introduces a linear operator S such that Sc(x) = f(x) , (3.23) where I define f(x) = £ m ( x ) , (3.25) so that Equations 3.22 and 3.23 are equivalent. If I can find a function G(x,x') that satisfies SG(x, x') = 6{x - x') , (3.26) then the desired solution can be expressed as /oo G{x,x')f(x')dx' . (3.27) •oo Here S(x — x') is the Dirac delta function and x' is a parameter. It is important to note that in Equation 3.27, Green's function is defined on an infinite domain (—oo, co). Since my mathematical problem is described on a finite domain [0, L], I either search for a different Green's function, or an additional assumption regarding the boundary conditions must be introduced in order to be able to use Green's function on an infinite domain. For example, as seen in Lee et al. (2001), if the domain's boundaries 36 have little or no effect on the chemical distribution, then I can reason that the desired finite domain Green's function is just a small perturbation of the function defined for the infinite domain (Morse and Feshbach 1953). In addition, in my biological system the boundary conditions represent an abstraction. Equation 3.26 states that G(x,x') satisfies SG(x,x') = 0 for all x except for x = x'. At x = x', I have an infinite spike concentration. Thus, I am left with the task of finding G(x, x'). I solve Equation 3.26 in the two regions x < x' and x > x' where SG(x,x') = 0. (3.28) A linear equation such as 3.28 has exponential solutions G(x,x') = e rx, where r is the eigenvalue determined from the characteristic equation - ( r 2 - A ) = 0 . (3.29) Equation 3.29 yields two roots r = so the general solution to 3.28 is of the form G{x, x') = Vxe^x + Vie~^x , for - oo < x < x' , (3.30a) and G{x, x') = U2e^x + V2e~^x , for x' < x < oo , (3.30b) where Vj and Vj are constants. To determine all the coefficients V j and Vj , I use several properties of the function G(x,x'). From the boundary conditions, in region 1 (—oo < x < x'), G(x,x') —>• 0 as x —> —oo. Therefore, in region 1, I must have V\ = 0. Similarly, in region 2 (x1 < x < oo), G(x, x') —>• 0 as x —> oo, which implies U2 — 0. So far I have G(x, x1) Vie V D> ; for - oo < x < x' l r£ (3-31) ~ \ D X V2e V ; ; for x < x < oo . 37 To find Ui and V2, I use two additional conditions that must be satisfied by Green's function. The first one is a matching condition at x', i.e. G(x, x') is continuous at the point x', so that L I^-x' -L / 6' x' JJxe V D - = V2e VD< . (3.32) The second condition to be satisfied by G(x, x') is the so called jump condition at x' and it can be derived from Equation 3.26. I have /oo poo SG(x, x') dx= ' 6{x - x1) dx . (3.33) oo J —oo If I integrate the left-hand side across the point x' and use /oo 5(x — x') dx — 1 , -oo I get /x px T poo SG{x,x')dx+ SG{x,x')dx+ SG(x, x') dx = 1 , (3.34) -oo J X ~ JX + where x ' - and x ' + indicate that I consider only values of x' that are less than or greater than x', respectively. I know that SG(x,x') = 0 for all x except for x = x', so that the first and third integrals on the left-hand side vanish; using Equation 3.24, I obtain ' + fd2G(x,x') b or dG(x,x') b_ D px T / G(x,x')dx =-1 . dx Since G(x, x') is continuous at the point x', the second term in the previous equation is also equal to zero, and I obtain the jump condition G'(x,x')\x=x>+ ~ G'(s,*')!«=*'- = - 1 • (3.35) 38 Using Equation 3.31, Equation 3.35 becomes _ V 2 ^ e - v J * ' - U ^ e ^ b * ' = - 1 . (3.36) Equations 3.32 and 3.36 form a system of two equations with two unknowns, U\ and V2, that can be easily solved, yielding D 1 . / A V2 = — F e V o 2 . M Thus, I find that the Green's function is given by X •K=e V^D 1 'e ^ y~o x , for — 00 < x < x' G(x,x') = { *VD ( 3 . 37 ) - L = e V / 5 " I ' e ~ v / s ' x , for x' < x < 00 , or, its equivalent form, G(x,x') = y j e - ^ x - x 'I. (3.38) For a given set of parameters, I show the graph of Green's function in Figure 3.1. Intu-itively, Green's function represents chemical concentration at the point x resulting from the production of this chemical by cells at the point source x'. Finally, using Equations 3.25, 3.27 and 3.38, the solution c(x) of Equation 3.22 follows poo dx) = -4= / e-V>-x'\m(x') dx' . (3.39) 2VWJ-00 This is the desired relationship between the cell distribution and the resulting chemical profile that the cells secrete. The form of the chemical Equation 3.39 indicates that at a point x, the chemical concentration c(x) is found by superimposing concentrations resulting from cells secreting at all points x', — co < x' < co. 39 Figure 3.1: Green's function for a set of parameters: x' = 1, b = 2, D = 8. 7 Numerical Simulations One of the mathematical tools that I use throughout the thesis to analyze the models is numerical simulation. Therefore, in this section, I describe the numerical method of discretization of the initial boundary value problem given by the partial differential equations 3.4, the boundary conditions 3.5, and the initial conditions 3.6, followed by computer generated numerical simulations. This particular numerical method was im-plemented to obtain the numerical simulations generated by X P P 1 and shown later in the section. 7.1 The Finite Difference Method for Initial Boundary Value Problems To solve the initial boundary value problem consisting of Equations 3.4, 3.5 and 3.6, I convert the partial differential equations to finite difference equations by replacing the derivatives with central difference quotients (Burden and Faires 1989). Due to the implementation technique of the resulting equations on X P P , I only have to replace the spatial derivatives in Equations 3.4 and 3.5, as the partial differential equations are solved XPP=Xterm Phase Plan software developed by Bard Ermentrout for dynamical systems. 40 in time by the "method of lines". The domain of the boundary value problem is the interval [0,L]. I used a uniformly spaced grid consisting of 51 points; it is conventional to let the letter h stand for the uniform difference in the x-values, h = Ax, and in particular, here h = L/50. In Appendix A3, I derive formulas for the central difference approximations for the first and second order derivatives that appear in the right-hand side of Equations 3.4. Using the common notation for numerical approximations, I obtain /; = fn+l~hfn-1+0{h2)1 (3.40a) where f'n = f'(xn) and f£ = f"(xn) are the generic first and second order derivatives that will replace the derivatives in my model equations. / „ = f(xn) is given by the initial conditions. The order relation 0(h2) signifies that the error is proportional to h? as h —> 0 (Burden and Faires 1989). Equations 3.40 are called central difference formulas for the first and second order derivatives. Central difference approximations to derivatives are more accurate than forward or backward approximations (0(h2) error versus 0(h) error). For comparison, I have the following 0(h) approximation for the first derivative: f ^ = f n + 1 ~ f n + 0 ( h ) . (3.41) Formula 3.41 is known as the forward difference formula if h > 0 and the backward difference formula if h < 0. In addition to the central difference formulas given by 3.40, I also used a second order approximation for Equation 3.4a (R. Russell, personal communication). Using this ap-proximation, I do not expand the second term of the right hand side. Instead, I use a 41 midpoint approximation, as follows: m-dx \ dx) { h _ + rrij) (d+i - Cj) - (rrii + m^i) (ct - C j _ i ) 2 / i 2 When the boundary conditions involve derivatives, the domain must be extended be-yond the interval, using fictitious points (Gerald and Wheatley 1989). As shown in Appendix A3, this leads to the following discretized cell equation: ^ ( 0 ) = 2 ^ - X ( m i + m f 1 " C ° ) , (3.43a) dm mi+i - 2mt + m^i {mi+x + m i ) ( c i + i - Cj) - (rrij + mi-i)(cj - C j _ i ) dt [ l ) ~ 1 1 h? X 2h? i = l , . . . , 4 9 , (3.43b) d m ( z c \ \ - o m 49 ~ m 50 ( m 4 9 + m 5 0 ) ( c 4 9 - c 5 0 ) "aTC50 ;" ^ x tf • ( j Similarly, the chemical equation can be expressed as follows: '• ^(0) = 2D^^+g(m0,c0) (3.44a) f (0 = ^ Q + 1 ~ 2 ^ + C'"1 + gim, Ci) , i = 1,..., 49 , (3.44b) ^ ( 5 0 ) = 2 D ^ ^ + s(m 5 0 , c 5 0) . (3.44c) The set of Equations 3.43 and 3.44 constitutes the full discretization needed to implement System 3.4 on X P P . A couple of important remarks must be made. First, smaller values for h should improve the precision. Second, I have to be concerned with stability and convergence criteria 42 (Gerald and Wheatley 1989). I define the following ratio (3.45) where D is the diffusion coefficient of the chemical, h = Ax is the space step size, and A t is the time step size chosen for solving System 3.4. In order to ensure stability and convergence in the explicit method, I must have which represents the Courant-Friedrichs-Lewy condition (Press et al. 1989). By conver-gence is meant that the numerical results of the method approach the true solution of the partial differential equation as A t and h both approach zero. Stability is related to propagation of errors. Thus, errors made at one step of the calculations should not cause increasingly large errors as the computations are continued, but rather will eventually damp out. 7.2 The Diffusion Equation In Section 5, I analytically calculated the chemical solution of the diffusion Equation 3.14, as given by 3.18. Here I use computer generated numerical simulations to validate the fact that, as t —» oo, the chemical concentration should approach its homogeneous steady state c, as determined by the initial conditions. In addition, I would like to investigate the influence of the diffusion coefficient D on the chemical distribution. Using the form of the chemical Equation 3.18, I propose the following initial condition: where n is the mode and L is the domain length arbitrarily set at L = 50. Clearly, the homogeneous steady state is c = 1. (3.46) x, 0) = 1 + cos(——) , (3.47) 43 0 I 1 1 1 1 1 0 500 1000 1500 2000 2500 time Figure 3.2: Time behaviour of the chemical concentration c(x2i,t) in the diffusion Equation 3.14, for different values of the diffusion coefficient D. The initial condition is c(x,0) = 1 + cos(47ra;/L). Simulations generated on XPP. Figure 3.2 shows the influence of the diffusion coefficient D on the dynamics of the Equation 3.14. I consider three different values for the parameter D, i.e. D = 0.05, D = 0.5, and D — 2, keeping the initial condition the same. I expect, and observe, that the higher the value of D, the faster the homogeneous steady state c = 1 is reached. 8 Conclusion In this chapter, I introduced the mathematical background needed for the model. I first defined the modelling goals and then presented the general system used to describe my chemotactic model, in an arbitrary coordinate system, followed by the one dimensional system. While analyzing subsets of the one dimensional chemotactic system, I also introduced several mathematical tools that will be used throughout the thesis. Thus, I showed how to find stationary solutions of a system; I presented the method of Green's function which can analytically solve partial differential equations, if certain conditions are satisfied by 44 the system; I generated numerical simulations. In general, numerical simulations are helpful for at least two reasons: first, they can corroborate theoretical results; second, they suggest new directions of study. 45 Chapter 4 The Keller-Segel Chemotactic Model In this chapter, I present a detailed analysis of the Keller and Segel (1970) chemotactic model introduced first for the slime mold aggregation. I chose to examine this particular model because it is the simplest mathematical model for aggregation of cells due to chemotaxis in response to autonomously created chemical concentration gradients, and therefore I expect to gain analytical insight from its investigation. My goal is to determine what types of interactions between cells and chemicals would lead to the formation of aggregates in the Alzheimer's disease brain. The model equations and the mathematical tools that I use to analyze the chemotactic model are shown in the first section. Several simplifications and reductions of the system follow. I describe in detail the linear stability analysis and the derivation of the sta-tionary solutions for this model, and I conclude with numerical simulations. Biological interpretations of the model are also discussed. 1 Model Equations and Mathematical Tools Recall that I model the interactions between microglia cells and toxic chemicals using the one dimensional chemotactic system of equations 3.4. I defined m(x, t) to be the density of activated microglia at time t and c(x, t) to be the concentration of the toxic chemical, where 0 < x < L and L is the domain length representing a typical dimension in the 46 brain. c(x, t) can be defined as the concentration of /3-amyloid protein or, alternately, of some chemical secreted by microglia. The neurons are not explicitly considered, but it is assumed that locally high levels of c(x, t) are representative of putative plaques, and are injurious to neurons. Microglia movement is characterized by two parameters (assumed constant), the random motility \i and the chemotaxis coefficient x- In this model, I consider that the cells do not proliferate or die. The chemical factor can diffuse in the region with constant diffusion coefficient D. The term of interest in the model is the chemical kinetics function, g(m,c) which is assumed to be the linear function g(m, c) = am — be. This reaction function indicates that microglia produce chemicals at a constant rate a, while the chemical decays at a constant rate b. I chose this form of g(m,c) because it is the simplest function that models the secretion of attractant by cells. I obtain the following system: dm d2m d ( dc\ , A , \ dc d2c « = Dw+am-bc- (4-lb) As introduced in the previous chapter, I assume that the problem has no flux boundary conditions, i.e. dm dx dc dx = 0 , (4.2a) x=0,L = 0 , (4.2b) x=0,L so that cells and chemicals do not leave or enter the domain. The initial conditions are m(x, 0) = u(x) , 0 < x < L , (4.3a) c(x, 0) = w(x) , 0 < x < L , (4.3b) 47 where u(x) and w(x) are some sufficiently smooth functions. The mathematical tools that I use to analyze this model are: • linear stability analysis that provides conditions for growth/decay of initial aggre-gates (Keller and Segel 1970). • analysis of the stationary system which proves existence of non-homogeneous sta-tionary distributions (Schaaf 1985; Schaaf 1984). • numerical simulations to identify stable non-homogeneous solutions. 2 Simplifications and Reductions In order to analyze System 4.1, the following two reduced versions of this system will be considered. Each of these subsystems will lead to certain properties, conditions and results to be further explored. 2.1 Space-Independent System Biologically, this is equivalent to looking at a well mixed, uniformly distributed system of microglia and chemical. Mathematically, a space-independent system is the system in which there are no spatial gradients, i.e. This system will give rise to spatially uniform steady states of the System 4.1. dm ~dt dc dt am — be . 0, (4.4b) (4.4a) 48 2 .2 Time-Independent System A time-independent system is a system in which there are no time derivatives. The solutions of such a system are spatially non-uniform steady states and they satisfy d2m d ( dc\ 0 = ^ - ^ r ^ j ' ( 4 - 5 a ) d2c 0 = D — + am-bc. (4.5b) ox2 We are already familiar with finding stationary solutions for a system, as introduced in Chapter 3, Section 4.1. 3 Linear Stability Analysis of the Space-Independent System In this section, I examine the stability of the homogeneous steady states of the System 4.1, using the linear System 4.4. The study is known as the linear stability analysis. My model must have non-trivial solutions in order to be biologically relevant. Thus, if I let the homogeneous steady state be (m, c), then it must satisfy m(x,t) = constant, (4.6a) am-be = 0 . (4.6b) For a given set of initial conditions, the value of m is determined uniquely (by conserva-tion) from m= m(x, 0) dx . (4.7) Jo The expression for c is then c = -m . (4.8) 49 Stability of the homogeneous steady state (m, c) depends on the Jacobian J of the System 4.4, ' 0 0 (4.9) \ a - b j I define T = trace J = -b , (4.10) A = d e t J = 0 . (4.11) Since A = 0 and r < 0 for any values of the parameters a and b, the homogeneous steady state (m, c) is a stable node (see Appendix B l ) . In addition, the theory developed by Keller and Segel (1970) states that instability of this homogeneous steady state can lead to the onset of aggregation in a system. This becomes the subject of my analysis in the next section. 4 Linear Stability Analysis of the Full System I am interested in further analyzing the possible change of stability of the homogeneous steady state (m,c). From Keller and Segel (1970) I know that spatially non-uniform perturbations of a homogeneous steady state could give rise to pattern formation. If the perturbations grow then I obtain peaks that represent aggregation of cells. In a neighbourhood of the homogeneous steady state (m,c), I let m(x,t) = m + m'(x,t) , (4T2a) c(x,t) = c + c'(x,t), (4.12b) where m'(x,t) and c'(x,t) are small perturbations from the steady state. I substitute 4.12 into System 4.1 to obtain dm' d2m' f dm'dc' . <92c'\ / „ , „ . V^rr - X + (™ + m , (4.13a) dt dx2 \ dx dx dx2dc' d2c' °^ = D ^ + am'-bc'. (4.13b) 50 Since the perturbations m! and d are small, the terms that are quadratic in the pertur-bations or their derivatives are extremely small, and I therefore neglect them to find the linearized system dm d2m _ d2c , A , A \ ~dl = ( 4 1 4 a ) dc d2c — = D— + am-bc , (4.14b) where I dropped the prime signs for simplicity. Equations 4.14 are linear second order parabolic partial differential equations. Therefore, they have eigenfunction solutions in the form of sines and cosines. I thus look for solutions of the form m(x,t) = Meat cos(qx) , (4.15a) c(x,t) = Ceatcos(qx) , (4.15b) where • q is the wavenumber. In Appendix A l , I have shown that q = rnr/L. The positive integer n is called the mode. Biologically, ir/q represents the distance between peaks, i.e. the spacing of the aggregates. • cr, the eigenvalue, represents the growth rate of the perturbations. Its sign will determine whether I'll obtain growth or decay of perturbations. • M and C are amplitudes of the perturbations at time t = 0. Substituting 4.15 into System 4.14, I obtain a system of algebraic equations in M and C: (a + iiq2)M -xfhq2C = 0 , (4.16a) -aM + (a + Dq2 + b) C = 0 . (4.16b) 51 System 4.16 has a trivial solution M = C = 0. Non-trivial solutions can be found if and only if the determinant of the System 4.16 equals zero, i.e. ( a + uq2 — ymg 2 \ = 0 . (4.17) -a a + Dq2 + b J It follows that (a + pq2)(a +Dq2+ b) - xmaq2 = 0 , (4.18) which leads to the d i s p e r s i o n r e l a t i o n given by a2 + 8a + J = 0, (4.19) where 8 = q2(fi + D) + b, (4.20a) 7 = q2 [fi(Dq2 + b) - Xma] . (4.20b) The quadratic Equation 4.19 has two roots: °i,2 = g ' ^ ' Searching for perturbations that would cause instability of the steady state, I look for conditions such that Real(a) > 0. Clearly 8 > 0 for any values of the parameters, and therefore one of the roots always has a negative Re(c). To obtain a second root with positive real part, I need 7 < 0. This is equivalent to Xfna > n(Dq2 + b) . (4.22) Inequality 4.22 represents the a g g r e g a t i o n c o n d i t i o n or the l i n e a r i n s t a b i l i t y c o n -d i t i o n . The inequality involves the wavenumber q and several other parameters. For example, decreasing /i, D or b, and/or increasing X o r a i i - e - having lower cell motil-ity, diffusion or degradation rate of chemical, or higher cell chemotactic sensitivity or 52 secretion rate of chemical, makes it more likely that Inequality 4.22 would be satisfied. Thus, whenever Inequality 4.22 holds, the homogeneous steady state (m, c) is unstable to perturbations, and aggregation of microglia occurs. Figure 4.1 shows the dispersion relation 4.19 for a given set of arbitrarily chosen param-eters. Aggregation of cells occurs in the region where a > 0. I observe that by i n c r e a s i n g a, the secretion rate of the chemical, the range of possible values of q also increases. This implies that the mode n could be bigger, and therefore, the number of aggregates per domain length should increase. (Actual results would depend on non-linear interactions not considered here.) q Figure 4.1: Dispersion relation for the Keller-Segel model, for different values of a. Parameters: b = 1, D = 1, x = 0-1, n = 1, rh = 10.8. 5 Biological Interpretations of the Aggregation Con-dition We have seen in the previous section that whenever the linear instability condition given by Inequality 4.22 is satisfied, pattern formation occurs in the Keller-Segel system. There-fore, I can use this information to interpret biologically the aggregation condition, make a connection between the mathematical model and Alzheimer's disease, and maybe suggest 53 therapy interventions for the disease. Pattern formation in the Keller-Segel model can be correlated with plaque formation in Alzheimer's disease, since both form due to clustering of cells and chemicals. Biologi-cally, I am interested in conditions that would slow or stop the formation of plaques in Alzheimer's disease. Mathematically, if the aggregation condition is N O T satisfied, ag-gregates will not form. Therefore, using Inequality 4.22, I look for conditions that could lead to NO plaque formation in Alzheimer's disease. I suggest the following conditions that could prevent or affect plaque formation in Alz-heimer's disease (one may speculate that these suggestions would be beneficial, but these ideas may be controversial since many Alzheimer's disease experts disagree on what are causes and what are mere downstream effects of the pathology): • decrease fh, the number of activated microglia. This research is already under investigation. McGeer and McGeer (1999) use anti-inflammatory drugs to reduce inflammation and decrease the proportion of activated microglia in the brain. • decrease a, the production rate of toxic chemical. A group of researchers at the Cen-ter of Neurological Diseases at Harvard's Brigham and Women's Hospital, headed by Denis Selkoe (Nash 2000), are studying this phenomenon. They test several drugs to inhibit the production of the /3-amyloid protein. (I mentioned at the be-ginning of this chapter that c(x, t) could represent /3-amyloid protein or any other chemical secreted by microglia.) • increase b, the removal rate of /3-amyloid protein. • decrease x> the chemotactic coefficient of microglia. • increase fj,, the microglia motility rate, and D, the /3-amyloid protein diffusion. 54 6 Stationary Solutions 6.1 The Cell Equation The form of the time independent System 4.5 allows for near-explicit stationary solutions. This has been described in an advanced mathematical treatment (Schaaf 1985) for the linear function g(m,c) = am —be. I use Schaaf's technique to understand several aspects of the system. The stationary system is given by Equations 4.5. Equation 4.5a can be rewritten in the form: dx dm dc dx dx (4.23) As shown in Chapter 3, Section 4.1, successive integrations of Equation 4.23, with the no flux boundary conditions 4.2, lead to the non-trivial solution m = m{c(x)) = M e * c ( x ) . (4.24) I note that this establishes a connection between microglia and chemical distribution which must hold whenever a non-uniform stationary solution of the original system exists. Furthermore, I observe that this relationship contains a free parameter, M. I can now use this result to eliminate the variable m from Equation 4.5b. 6.2 The Chemical Equation: The Reduced System With the result of the previous section, I can now rewrite Equation 4.5b in the form ol2c x 0 = D— + aMeic - be , (4.25) where the function to be identified is c(x). I observe that this equation is not linear in c, and therefore, I cannot solve it using the method of Green's function introduced in Chapter 3, Section 6.1. 55 To understand qualitatively how solutions of this second order ordinary differential equation behave, I define w = dc/dx and set up an equivalent system of two first or-der equations dc ,, n n . — = w , 4.26a dx d£= -iG{c), (4.26b) where I have defined the new function G(c) = a M e * c - be . (4.27) I refer to Equations 4.26 as the Reduced System. The function G(c) contains the arbitrary constant M , which will serve as a bifurcation parameter (Schaaf 1985). A l l non-uniform stationary solutions with non-trivial microglia cell density given by Equation 4.24 satisfy this system of equations. I can study the behaviour of this system in the cw phase plane and use the qualitative properties of the trajectories in that phase plane to understand the profiles of chemical concentration c, that can occur as stationary solutions. I also used the information generated by the Reduced System to implement my numerical simulations. 6 .3 Fixed Points of the Reduced System I would like to find the fixed points (c*,w*) of the Reduced System. Note that fixed points of the Reduced System represent homogeneous steady states of the original System 4.1, since the spatial gradient of the chemical c is set to zero in System 4.26. Thus, the fixed points of System 4.26 must occur at w* = 0 by Equation 4.26a. They must also satisfy G(c*) = 0 which is equivalent to (4.28) As shown in Appendix B2, the analytical solution of the transcendental Equation 4.28 can be derived using the method of Lamber t W function (Corless et al. 1996). The 56 solution c* is then given by c* = -^W(v) , (4.29) X where W(v) is the Lambert W function and v — — (XaM)/(fib). If — 1/e < v < 0, then there are two possible real values for W(v) that can be calculated using M a p l e (Corless et al. 1996). The arbitrary constant M is here treated as a bifurcation parameter. Thus, for the Keller-Segel chemotaxis model, there are at most two fixed points in the phase plane. The character of the fixed points organizes the trajectories in the cw plane, and therefore determines the profiles of m(x) and c(x) that could occur in stationary solutions. Stability properties of these fixed points play an important role. The Jacobian matrix associated with System 4.26 is J=\ , ° , 1 ) , (4.30) as follows from Appendix B l . I also calculate r = 0 , (4.31) A = 1 (jc- - l ) . (4.32) The trace r is zero, and therefore the fixed points of the Reduced System can only be classified as either saddle point or neutral stable center, depending on the sign of the Jacobian's determinant A . If A > 0, i.e. c* > p,/x, then the fixed point (c*,0) is a center. For A < 0, i.e. c* < pjx-, the fixed point (c*, 0) is a saddle. 6.4 The Phase Plane cw In Figure 4.2, I plotted the phase plane cw for a given set of parameters. Before I can use the qualitative properties of the trajectories in the phase plane cw in order to understand the profiles of chemical concentration c that can occur as stationary solutions of the Keler-Segel chemotaxis model, I note the following important points: 57 Figure 4.2: The phase plane for the Keller-Segel model. Parameters: a = 1, b = 1, Af = 3.6143, JD = 0.3, x = 0.1,/i = 1. 1. Trajectories in this cu; plane have as the independent variable the spatial variable x. Thus, I have to interpret the trajectories in terms of spatial profiles of chemical concentration over the domain. 2. To be biologically relevant, these trajectories must stay in the positive c semi-plane, and must be bounded. (The value of w, which corresponds to gradients of c can be negative, but it too, should be bounded for realistic solutions.) 3. To satisfy the boundary conditions given by Equations 4.2, only trajectories satis-fying w — 0 at x — 0, L are candidates, i.e. trajectories must start and end on the line w = 0. 4. Fixed points of this system, i.e. values which satisfy dc/dx = dw/dx — 0, represent profiles that are spatially homogeneous. In order to satisfy the constraints listed here, I must restrict my attention to the following types of trajectories, whose interpretation is as stated: Homoc l in i c trajectories These are trajectories that start and end at some fixed point (which must therefore be a saddle point). If such a trajectory is based at c = 0, 58 it could represent a single "peak" of protein inside the domain tapering off to low levels at the boundaries. In principle, this trajectory actually describes the behaviour of c(x) over the infinite domain — oo < x < oo, so applying this to a finite domain is an approximation. Periodic trajectories These are bounded trajectories that represent some periodic arrangement of chemical peaks. Thus, they are of particular interest in my problem, especially when they also satisfy the boundary conditions 4.2. Fixed points These would represent solutions whose chemical profile is constant, i.e. spatially homogeneous. I know that such solutions can be the basis of pattern formation provided they satisfy the conditions for instability given in the linear stability analysis of the full system. Thus, bifurcations from these solutions are of interest to us, in particular bifurcations from the neutrally stable center point. (I am not concerned with bifurcations from the saddle point since such a point is, by definition, unstable.) In Figure 4.2, there are two fixed points: (c^,0) is a saddle point, while (c^O) is a center. Here, each bounded trajectory satisfying c > 0 represents a possible stationary solution of the system. Also, one of the fixed points (c*,0) and (c^O) represents the homogeneous steady state of Equations 4.1 (the homogeneous steady state is unique, as explained in Section 3); the circular orbits represent periodic patterns of alternating low/high concentration (as would occur in plaques) and the homoclinic orbit corresponds to a single cluster. 7 Numerical Simulations The numerical simulations shown in this section were generated for the Keller-Segel chemotactic model given by Equations 4.1, the boundary conditions 4.2 and the initial 59 conditions 4.3, using arbitrarily chosen parameters. They reproduce the cell and chemical dynamics predicted by the linear stability analysis done for the model. The simulations were obtained by using an adaptive mesh program developed by R. Russell, Simon Fraser University, and in collaboration with A. Chavez-Ross. The program can follow sharp spikes when the concentration of cells or chemical increases considerably within a given aggregate (Figures 4.6 and 4.7c,d,e). Also, it shows how periodic distributions decay (Figure 4.5a), or grow (Figure 4.6a), depending on the parameter values. 7.1 Derivation of Parameters In order to generate numerical simulations, I have to develop a step-by-step scheme to find sets of parameters that would (or would not) satisfy the aggregation condition 4.22. The scheme is as follows: Step 1 System 4.1 contains five parameters. I first choose the following arbitrary values for these parameters: a = l, 6 = 1, r> = 0.3, x = 0.1, /i = l . Step 2 I must also find m and c to complete the specification of the system. The idea is to select m and c such that for some integer n > 0, the initial conditions represent spatially non-uniform perturbations of the homogeneous steady state (m,c), perturbations that give rise to pattern formation. I choose e = 0.1 and the domain length is arbitrarily chosen as L = 50. The homogeneous steady state (m, c) and the mode n will be determined in the following steps. m(x, 0) (4.33b) (4.33a) 60 Step 3 I mentioned that bifurcations from the center fixed point of the Reduced Sys-tem 4.26 are of interest to us. I also determined that a center occurs when c* > /i/x-Therefore, I restrict my attention to values of the homogeneous steady state c such that c = c* > 10 for the above parameters. Step 4 The fixed point c* is given by Equation 4.29 and its numerical value can be exactly determined using Maple, if the value of the bifurcation parameter M is also provided. If I impose — 1/e < v < 0, it follows that xaM/[ib < 1/e. Given the parameters values in Step 1, I obtain M < 3.676. Let M = 3.6. Step 5 Using Equation 4.29, I find c* = 8.2 and c* = 12. Let c = 12 > 10. Step 6 I find the homogeneous level of microglia density m to be also fa — VI, either from Equation 4.8 or 4.24. So the homogeneous steady state of the original System 4.1 is set to be (m, c) = (12,12). Step 7 I rewrite the aggregation condition 4.22 in the form This inequality can be understood graphically, as shown in Figure 4.3. The left-hand side y — c — n/x describes a linear function of c with slope 1 and intersecting the horizontal axis at c = fi/x- The right-hand side is a constant for a fixed setting of the parameters, and a choice of the mode number n. For the given parameter values, the inequality is satisfied, but only for waves of mode n < 12. The unstable and stable modes n can also be visualized if I plot the dispersion Equation 4.20 (Figure 4.4). I have unstable modes when a > 0, and the mode n is calculated from q = n i r / L . The graph indicates that the most unstable mode n occurs at a m a x , i.e. the leading positive eigenvalue associated with the system of Equations 4.1. In this particular case, the most unstable mode could be either n = 8 or n = 9. (4.34) 61 3 2.5 2 0.5 0 / / 6 - \U% I / n = 14 / n = 13 J. 1 I n = 12 / n = 11 / _ / n = 10 / / / / / / C = 12 Figure 4.3: Plot of the linear function c — u/x and of the constants on the right-hand side of Inequality 4.34 for the Keller-Segel model. Parameters: b=l,D = 0.3, x = 0.1, M = 1,L = 50. q c = 0.535 1 Pmax = 0.815 0.4 0.6 q Figure 4.4: Plot of the dispersion Equation 4.20 for the Keller-Segel model. Pa-rameters: a = 6 = // = 1, D = 0.3, x = 0.1, m = 12, L = 50. The most unstable wave mode n occurs at qc. 62 Step 8 I also provide random initial conditions: m(x, 0) = m + ex ran > (4.35a) c(x, 0) = c + ex ran j (4.35b) where e = 0.002 and x r a n e [0,1] is a random number internally generated by the com-puter. A l l the other parameter values are as previously described. With these initial conditions it is difficult to make any predictions about the number of peaks forming in the domain since any mode n < 12 can grow and give rise to pattern formation. More importantly, the nonlinear terms that I neglected in my linear stability analysis become important when stationary solutions form in such a chemotactic system (Segel 1984). 7.2 Graphs Having chosen all the parameter values and the initial conditions, I produce two sets of numerical simulations: one set of graphs were obtained using periodic initial con-ditions (see Figures 4.5 and 4.6), while a second set of graphs were generated using random initial conditions (see Figure 4.7) Thus, for n = 14, the aggregation condition 4.34 is not satisfied and the system ap-proaches its homogeneous steady state, as can be seen in Figure 4.5. Figure 4.5a shows the decay of microglia distribution. In Figure 4.5b, I show the microglia spatial distri-bution at time t = 0. The graph reproduces the total number of 7 peaks predicted by using periodic initial conditions for mode n = 14. In contrast, for mode n = 10, instability occurs and stationary pattern of microglia develops, as shown in Figure 4.6a. In Figure 4.6b, I plotted both the microglia and the chemical spatial distributions at time t = 500. The number of 5 peaks validates the use of the periodic initial conditions with wave of mode n — 10. In Figure 4.6b, we can also observe the exponential relationship established between microglia and the chemical, as 63 (a) Decay of the microglia distribution. (b) Microglia distribution at time t = 0. Figure 4.5: Microglia distribution in the Keller-Segel model initiated with the sta-ble mode n = 14. Bifurcation parameter M = 3.6143. Parame-ters: a = b = u = 1,D = 0.3, x = 0.1, e = 0.1, L = 50. Ho-mogeneous steady state: rh = c — 12. Periodic initial conditions: m(x, 0) = m + e cos(nir/L), c(x, 0) = c + e c o s ( n n / L ) . predicted by the existence of stationary solutions of the Keller-Segel chemotactic system. The second set of graphs were obtained using random initial conditions. Figure 4.7 presents the time evolution of both microglia and chemical distributions, from time t = 0 to time t — 500. At time t = 0, I see in Figure 4.7a the random spatial distribution of the cells. Later, at time t = 200 (see Figure 4.7b), 4 peaks have evolved. At time t = 420 (see Figure 4.7c), 3 of the peaks have decreased to very low relative sizes and there is only one big peak of microglia forming. Finally, at time t = 500 (see Figure 4.7d), 4 sharp peaks of microglia and 4 lower peaks of chemical form. The number of peaks correlates with the most unstable mode n = 8 that I found by plotting the dispersion relation in Figure 4.4. Several simulations were generated to ensure accuracy; they led either to Figure 4.7d or to Figure 4.7e. Figure 4.7e shows formation of 4.5 peaks; this 64 (a) Decay of the microglia distribution. 12.15 | , , 1 , 11.85 ' 1 1 1 ' 1 0 10 20 30 40 50 space (b) Microglia distribution at time t = 0. Figure 4.5: Microglia distribution in the Keller-Segel model initiated with the sta-ble mode n — 14. Bifurcation parameter M = 3.6143. Parame-ters: a = b = u = l,D = 0.3, x = 0.1, e = 0.1, L = 50. Ho-mogeneous steady state: m = c = 12. Periodic initial conditions: m(x, 0) = m + e cos(n7r /L), c(x, 0) = c + e cos(nir/L). 64a •• microglia \ A A' j , chemical / \ \ / ' IX i 1 ; ! W-\ \ • V VI \1 \i \j • i 200 (a) Microglia and chemical (b) Microglia and chemical distributions at time t = 0. distributions at time t = 500. Figure 4.6: Stationary pattern formation in the Keller-Segel model initiated with the unstable mode n = 10. Bifurcation parameter M = 3.6143. Pa-rameters: a = b = /J, = 1,D = 0.3, x = 0.1, e = 0.1, L = 50. Ho-mogeneous steady state: m = c = 12. Periodic initial conditions: m(x,0) = m + ecos(n7r/L),c(a;,0) = c + ecos(mr/L). corresponds to the most unstable wave of mode n = 9. Thus, my numerical results clearly suggest that spatially non-uniform random perturbations of the homogeneous steady state of the Keller-Segel chemotactic system give rise to pattern formation with wave of mode n, where n is the most unstable mode calculated for a given set of parameters. 8 Conclusion In this chapter, I analyzed in detail the one dimensional Keller-Segel chemotactic system that models cell movement in response to chemical gradients. I investigated linear sta-bility analysis; I studied the stationary solutions of the system using several techniques applicable in the phase plane; I concluded with numerical simulations. As mentioned at the beginning of the chapter, by studying the Keller-Segel chemotactic model, I gained analytical insight of the system; this will become obvious in the next chapters when I will often apply similar analyses to new models. 65 space (a) Microglia and chemical distributions at time t = 0. 450 400 350 300 i 250 *, 200 E 100 50 IV. microglia chemical / space (b) Microglia and chemical distributions at time t = 500. Figure 4.6: Stationary pattern formation in the Keller-Segel model initiated with the unstable mode n = 10. Bifurcation parameter M = 3.6143. Pa-rameters: a = b = /j, = l,D = 0.3, x — 0.1, e = 0.1, L = 50. Ho-mogeneous steady state: rh = c — 12. Periodic initial conditions: m(x, 0) = rh + e cos(n7r /L) , c(x, 0) = c + e cos(mr/L). 65a (c) t = 420. (d) t = 500. microglia chemical V . . . h / ' \. J. 0 10 20 30 40 60 (e) t = 500. Figure 4.7: Spatial distributions of microglia and chemical in the Keller-Segel model, using random initial conditions. Bifurcation parameter M = 3.6143. Parameters: a = b = n = l,D = 0.3,% = 0.1,e = 0.002,L = 50. Homogeneous steady state: fh = c = 12. Random initial condi-tions: m(x, 0) = m + ex r a n , c(x, 0) = c + exT&n. 66 Numerical simulations using periodic initial conditions complement the linear stability analysis. Small spatially non-uniform perturbations of the homogeneous steady state could either decay or lead to the formation of stationary periodic aggregates of cells and chemical. In general, having lower cell motility, diffusion or degradation rate of chemical, or higher cell chemotactic sensitivity or secretion rate of chemical, makes it more likely that aggregates would develop in the chemotactic system. In addition, my numerical simulations using spatially non-uniform random perturbations of the homogeneous steady state indicate that the Keller-Segel chemotactic system shows periodic pattern formation with the most unstable mode n associated with the leading positive eigenvalue of the system. I note that the linear kinetics function g(m, c) = am — be is probably too simplistic to reproduce real biological interactions that exist between the cells in the brain. This function models a linear rate of secretion of chemical by microglia, and does not represent the removal or uptake of certain chemicals by these cells, nor any saturation phenomena. Therefore, in the next chapter, I will attempt to extend my model by considering different chemical kinetics function, and also by introducing a kinetics function in the cell equation. 67 Chapter 5 Extensions of the Keller-Segel Model In this chapter, I extend my Keller-Segel chemotactic system of equations to include other kinetics functions in both equations. There are two types of extensions that I study. One derives from the desire to make the mathematical models biologically meaningful. I accomplish this by considering a more relevant kinetics function in the chemical equation. I also incorporate a kinetics function in the microglia equation and investigate its effects on the dispersion relation. 1 Different Kinetics Functions for the Chemical Equation The one dimensional system of equations that models the interactions between microglia and chemicals introduced in Chapter 3 is: (5.1a) d2c + g(m,c) . (5.1b) dx2 I assume that the problem has no flux boundary conditions, i.e. dm dx x=0,L = 0 (5.2a) dc dx x=0,L = 0 . (5.2b) 68 The initial conditions are m(x, 0) = u(x) , c(x, 0) = w(x) , 0 < x < L , 0 < x < L . (5.3a) (5.3b) The kinetics function g(m,c) in Equation 5.1 plays a crucial role in the dynamics of the system. The study of this chemotactic system led to the investigation of several forms of the function g(m,c) which would be biologically relevant for studying senile plaque formation in Alzheimer's disease. I first analyzed the linear reaction function g(m,c) = am — 6c, as seen in the previous chapter. I would like to study more general forms, such as g(m, c) = s(c) — mh(c) , (5.4) where s(c) is the microglia-independent rate of chemical production in the tissue (e.g. by neurons), and h(c) is a function that represents the rate at which microglia either produce or remove the chemical, or both. I also require that chemical production terms are self-limiting, i.e. that s(c) and h(c) decrease as c increases. a. The function s(c) = ac(l — c). b. The function h(c) = bc(cCT\t — c)e Figure 5.1: Sketches of the functions s(c) and h(c) in g(m,c) = s(c) - mh(c). Parameters: a = 1, b = 1, ccrit = 2.5, a = 1. Thus, I proposed s(c) = ac(l — c) (see Figure 5.1a), and I considered two study cases: 69 0.4 I 1 i 1 1 1 r 0.3 r--0.4 I 1 1 1 1 1 1— 0 0.2 0.4 0.6 0.8 1 1.2 C a. The function s(c) = ac(l — c). 0.7 b. The function h(c) = 6c(cc ri t — c)e a c . Figure 5.1: Sketches of the functions s(c) and h(c) in g(m,c) = Parameters: a = 1, b = 1, c c r j t = 2.5, a = 1. 69a Case 1 h(c) = b(l — c) represents a linear production rate of chemical. For this kinetics function, the level c = 1 corresponds both to the steady state chemical concen-tration in absence of microglia and to the level at which microglia switch from being removers to being producers of chemicals. I also observe that if c is very low, but positive, the cells remove the chemical at the roughly constant rate b per cell. Furthermore, in regions where microglia density is very large and c > 1 there could be an unlimited build-up of chemical. I studied this function extensively and the results were very similar to the ones already shown for the linear function g(m, c) = am — be. Therefore, this case will be omitted from further analysis. Case 2 h(c) — 6c(c c r i t — c)e~ac (see Figure 5.1b). Here h(c) represents a limiting pro-duction rate of chemical, c = c c rj t is the level at which microglia revert from being removers to being producers of the chemical, b is the rate of microglia up-take/absorption, and a represents the rate at which microglial secretion shuts off at high enough chemical concentration. (This shut off mechanism is required to prevent explosive concentration build-up at sites of aggregation of microglia.) In this case, g(m, c) = ac(l — c) — bc(ccrit — c)e~ac is the l i m i t i n g kinetics function. In addition, c = 1 is the level at which the production of the chemical factor turns off, and a is the rate of chemical production. In the next section, I introduce the conditions that led to the construction of the limiting kinetics function g(m, c). A detailed analysis of how the general features of this function influences the dynamics of my system then follows. 1.1 Conditions Imposed on the Kinetics Function g(m, c) The analysis of the Keller-Segel chemotactic model showed the existence and formation of stationary aggregates of microglia with very high densities. To account for this limitation, and additional biological implications, I assumed that the functions s(c) and h(c) have 70 the following features: 1. s(0) = h(0) = 0. This means that when there is no toxic chemical, the neurons are in good health and no toxin is being produced. Similarly, microglia are not removing the substance when it is absent. (This condition is necessary to avoid kinetics in which negative concentrations could occur.) 2. For small c > 0, s(c) > 0 and h(c) > 0. This means that the presence of some of the chemical factor, c, induces more production by the tissue and that microglia are acting as removers. 3. There is some steady state chemical concentration in the absence of microglia, and I call this level c = 1 without loss of generality. (This is equivalent to assuming that s(l) = 0 and s(c) < 0 for c > 1.) 4. There is some level cc r;t at which h(c) changes sign. For values of the chemical concentration c above this critical level, the microglia produce more toxic factor(s) than they consume, contributing to an increase in the local toxicity. 5. As the chemical concentration increases, h(c) should approach zero in order to level off production of toxic chemical. This is accomplished by using the function e~ac in the expression of h(c). 1.2 Dynamics of the Chemotactic System I now explore how the general features of the limiting kinetics function g(m, c) = ac(l — c) — bmc(ccv\t — c)e~ac influence the dynamics of the System 5.1. Thus, I ask: what parameter values give rise to pattern formation? I first identify the aggregation inequality associated with the System 5.1. I let the ho-mogeneous steady state (rh, c) satisfy g{m,c) = 0, (5.5) 71 or, its equivalent form, m He) (5.6) provided h(c) / 0. The linear stability analysis of the space-free system associated with the original System 5.1 shows that, in order to have a stable homogeneous steady state (m, c), I must have gc < 0, where gc = dg/dc(m, c) (see Appendix Cl .2) . Then, as shown in Appendix C l . 3 , the linear stability analysis of the full System 5.1 leads to the following dispersion relation: a2 + pa + 7 = 0 , (5.7) where ft = fiq2 + D q 2 - g c , 7 = <l2{^(Dq2 - gc) - xm% (5.8) (5.9) Here q — nir/L is the wavenumber and gm = dg/dm(m,c). Instability of the homoge-neous steady state (m, c) occurs when a > 0 which leads to the aggregation condition, or the linear instability condition: Xmgm > fi{Dq2 - gc) . (5.10) It is known that gc < 0, and therefore gm > 0 which is equivalent to c — c c rit > 0. Thus, for pattern to occur, I must have c c r i t < c in addition to the aggregation Inequality 5.10. Inequality 5.10 can be rewritten explicitly in terms of the functions s(c) and h(c) as follows (see Appendix Cl .3) : s c - s c h'(c) | X h{c) n\ >Dq2 (5.11) or its equivalent form ac 1 ccrit . /-. -\ ( X = + (1 - c) la Ccrit ~ c \ A4 > Dq2 . (5.12) 72 It is evident that the aggregation Inequality 5.12 can be satisfied provided that certain relationships between the parameters involved are satisfied. Recall that inequalities c c r i t < c and gc < 0 are also required to be satisfied. In the next section, I present a series of numerical simulations, generated for a particular set of parameter values, illustrating the dynamics of the chemotactic system with the limiting kinetics function g(m,c) = ac(l — c) — bmc(ccrit — c)e~ac. 1.3 Numerical Simulations The numerical simulations shown in this section were generated for the extended Keller-Segel chemotactic model given by Equations 5.1, the boundary conditions 5.2 and the initial conditions 5.3, using arbitrarily chosen parameters. I choose the following values for the parameters: and I let c = 3 > c c r j t . The microglia homogeneous steady state is calculated from Equation 5.6 to be m = 3.585. I verified that gc < 0, in particular gc = —3.6. I also provide periodic initial conditions: where the domain length is arbitrarily chosen as L = 20. e will take different values, depending on each particular simulations. These initial conditions represent spatially non-uniform perturbations of the homogeneous steady state (rn,c), perturbations that give rise to pattern formation, if additional conditions are also satisfied. The mode n in Equations 5.14 is determined from Figure 5.2 where I plotted the positive eigenvalue a calculated using Equation 5.7. Thus, for n < 9, the small perturbations Oi = Ccrit = 0.5 (5.13) (5.14b) (5.14a) 73 0.3 -0.05 I ' ' 1 1 ' ' ' '— 0 0.2 04 0.6 0.8 1 t.2 1.4 1,6 q Figure 5.2: Dispersion relation for the extended System 5 . 1 , with the limiting ki-netics function. Parameters: a = b = D = x = tJ' = ^ i a = c crit = 0.5. should increase in amplitude, creating instabilities in the system and giving rise to pat-tern formation. In contrast, when n > 10, the perturbed system should approach its homogeneous steady state ( m , c). I note that the most unstable mode is n = 6. I also used random initial conditions: m(x,0) = fh + exT&n, (5.15a) c(x,0) = c + e x r a n , (5.15b) where x r a n e [0,1] is a random number internally generated by the computer. A l l the other parameter values are as previously described. Thus, for n = 10, the aggregation condition 5.11 is not satisfied and the system ap-proaches its homogeneous steady state, as can be seen in Figure 5.3. Figure 5.3a shows the decay of microglia distribution. In Figure 5.3b, I show the microglia spatial distri-bution at time t = 0. The graph reproduces the total number of peaks (5) predicted by using periodic initial conditions for mode n = 10. In contrast, for mode n = 9, instability occurs and a stationary pattern of microglia (and chemical) develops, as shown in Figure 5.4a. In Figure 5.4b-f, I plotted the microglia 74 0.25 0.2 0.15 0.1 0.05 0 -0.05 1 1 • i i i 1 1 X — -y% = 0 - 9 8 ^ \ : y Pmax = 1-55 \^ i i i i i i i t 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 q Figure 5.2: Dispersion relation for the extended System 5.1, with the limiting ki-netics function. Parameters: a = 6 = D = x = /J = l ,a: = c c r i t = 0.5. 74a (a) Decay of the microglia distribution. (b) Microglia distribution at time t = 0. Figure 5.3: Microglia distribution initiated with the stable mode n = 10. Param-eters: a = 6 = D = x = / i = l , a = c c r i t = 0.5, e = 0.0001, L = 20. Homogeneous steady state: (m, c) = (3.585,3). Periodic initial condi-tions: m(x,0) = m + ecos(nn/L),c(x,Q) = c + e cos(nir / L). spatial distributions at different times. The time evolution of microglia density is different than the one obtained in the Keller-Segel model. At time t = 0 (Figure 5.4b), microglia forms 4.5 peaks, as expected, since I used periodic initial conditions with wave of mode n = 9. The number of peaks can still be observed at time t = 30 (Figure 5.4c), but at time t = 45 (Figure 5.4d), there are only 3 peaks which correspond to the most unstable mode n = 6. Much later, at time t = 480 and t = 1500 (Figure 5.4e and f), only one peak of microglia forms, a peak that sharply increases in density between the two given times. While generating my simulations, I found that for mode n = 10, instability can develop if e is not small enough. Thus, for e = 0.001, I obtained a similar spatial evolution of microglia density as for the unstable mode n = 9, as shown in Figure 5.5. In addition, when I used a smaller chemical concentration, i.e. c = 2, the system was not as sensitive to the value of e as in the previous case. For c = 2, even e = 0.01 generated the expected 75 1500 T space 3.5901 3.59008 3.59006 3.59004 3.59002 3 59 3.58998 3.58996 3.58994 3.58992 3.5899 3.58988 (a) Decay of the microglia distribution. V \J space (b) Microglia distribution at time t = 0. Figure 5.3: Microglia distribution initiated with the stable mode n = 10. Param-eters: a = 6 = D = x = M = l ) « = C c r i t = 0.5,e = 0.0001, L = 20. Homogeneous steady state: (m, c) = (3.585,3). Periodic initial condi-tions: m(x, 0) = m + ecos(n7r/L), c(x, 0) = c + e cos(nir/L). 75a (a) Microglia pattern formation. (b) t = 0. / \ A A ' \ / \ j \ / \ / \" v / \ / • \ / (c) * = 30. (d) t = 45. (e) t = 480. (f) t = 1500. Figure 5.4: Microglia distribution initiated with the unstable mode n = 9. Param-eters: a = b = D = x = V = l,a = c c r it = 0.5, e = 0.0001, L = 20. Homogeneous steady state: (m, c) = (3.585,3). Periodic initial condi-tions: 771(3:, 0) = m + e cos(n7r/L), C(X, 0) = c + e cos(mr/L). 7G microglia (a) Microglia pattern formation. (b) t = 0. (c) t = 90. (d) t = 240. (e) t = 2370. (f) t = 3000. Figure 5.5: Microglia distribution initiated with the "stable" mode n = 10. Pa-rameters: a = b = D = x = H = l,a = ccrit = 0.5, e = 0.001, L = 20. Homogeneous steady state: (m,c) = (3.585,3). Periodic initial condi-tions: m(x,0) = m + ecos(mr/L),c(x,0) = c + e cos(mr / L ) . 77 microglia 250 r 200 \-(a) Microglia pattern formation. (b) Microglia distribution at time t = 3500. chemical 9 . io r , . n space (c) Chemical pattern formation. (d) Chemical distribution at time t = 3500. Figure 5.6: Spatial distributions of microglia and chemical using random initial conditions. Parameters: a = b = D = x — a = 1, a = c c r i t = 0.5, e = 0.0002, L = 20. Homogeneous steady state: (m, c) = (3.585,3). Ran-dom initial conditions: m(x, 0) = fh + ea; ran, c(x, 0) = c + ex r a n . simulations. I conclude that the extended Keller-Segel model, given by Equations 5.1 with the limiting kinetics function, is very sensitive with respect to the amplitude of the small perturbations shown in Equations 5.14. The graphs shown in Figure 5.6 were obtained using the random initial conditions 5.15. 78 microglia 250 r 200 I (a) Microglia pattern formation. space (b) Microglia distribution at time t = 3500. 78a chemical 10 r space (c) Chemical pattern formation. 10 space (d) Chemical distribution at time t = 3500. gure 5.6: Microglia distribution initiated with the stable mode n = 10. Param-eters: a = & = D = x = M = l , a = Ccrit = 0.5, e = 0.0001,1 = 20. Homogeneous steady state: (fh,c) = (3.585,3). Periodic initial condi-tions: m(x, 0) = rh + e cos(n7r /L), C(X, 0) = c + e cos(mr/L). 78b The linear stability analysis does not predict how the chemotactic system evolves in this case, but my results are similar to the ones obtained for the unstable mode n = 9, so that the microglia and chemical distributions go through the same pattern changes as shown in Figure 5.4b-f. Thus, Figures 5.6a and c present the time evolution of both microglia and chemical distributions. At final time t = 3500 (see Figure 4.7c and d), one peak of each microglia and chemical form. I can also observe the exponential relationship established between microglia and the chemical, as predicted by the existence of stationary solutions of the extended Keller-Segel chemotactic system. 2 Multiple Kinetics Functions for the Chemotactic Model I propose a system of partial differential equations that includes both chemical and microglia kinetics functions g(m, c) and f(m, c), respectively. The function f(m, c) should incorporate cell division and/or death of microglia. The equations of this model are: din d2vn d ( dc \ •w = " & r - x s r s j + / ( m ' c ) ' ( 5- i 6 a ) | = D g + , < m , c ) , (5.16b) coupled with the boundary conditions 5.2 and the initial conditions 5.3. Systems like these were extensively studied by Maini et al. (1991) and Myerscough et al. (1998), and are known to produce periodic patterns. As I mentioned at the beginning of the chapter, this study was motivated by my search for a different dispersion relation than the one obtained in the Keller-Segel model. More specifically, I am investigating whether the eigenvalue a becomes positive for values of the modes n far away from zero. This scenario would lead to the formation of periodic patterns rather than large clusters (Cross and Hohenberg 1993). I summarize now the steps leading to the dispersion relation associated with my new 79 system of Equations 5.16 (detailed calculations are given in Appendix C2): • I introduce two reduced systems, the space independent system and the time inde-pendent system. • stationary solutions of the space free system represent homogeneous steady states (m, c) of the original system. • linear stability analysis leads to the dispersion relation and the aggregation condi-tion. Thus, the resulting dispersion relation is given by a2 + B a + 7 = 0 , (5.17) where B = w2 + D q 2 - f m - g c , (5.18) 7 = (iDq* - (xmgm + Dfm + iigc) q2 + fmgc-fc9m • (5.19) Here fm = df/dm(m,c), fc = df/dc(m,c), gm = dg/dm(m,c), and gc = dg/dc(m,c). Again, for instability to occur, I look for conditions such that Real(cr) > 0. Because there are too many parameters involved in the dispersion relation 5.17, I directly investigate it for different forms of the kinetics functions f(m, c) and g(m, c). Case 1 f(p,m) is logistic, g(m,c) is linear Let the kinetics functions be f(m,c) = r m ( l - , (5.20a) g(m,c) = am — be , (5.20b) where r is the microglia growth rate, k is the carrying capacity, a is the production rate 80 of the chemical, while b is the decaying rate of the chemical. It follows that B = (fx + D)q2 + b + r , 7 = pDq4 — (akx — rD — b/i)q2 + br . (5.21a) (5.21b) 0.3 0 0.2 0.4 0.6 0.8 1 1.2 1.4 q Figure 5.7: Dispersion relation for study Case 1 for several values of the bifurcation parameter b. Other parameters: a = 1, k = 1, r = 0.1, x = 3, p = 1, D = 1. The dispersion relation associated with System 5.16, for the kinetics Functions 5.20, is shown in Figure 5.7. The graph of a may or may not intersect the g-axis depending, in this particular case, on the bifurcation parameter b (all the other parameters are fixed). Since I am only interested in q > 0, the bifurcation occurs when a intersects the g-axis at some critical point, say qc, for a corresponding bc. Close to bifurcation, the linear stability theory then predicts periodic patterns of wavenumber qc (Segel 1984). I observe that by decreasing the decay rate of the chemical b, I can increase the range of unstable wavenumbers, though the linear stability analysis fails to predict which of these would eventually form the final (stable) pattern. 81 81a Case 2 f(p,m) is logistic, g(m, c) is limiting Let the kinetics functions be f(m,c) = r m ( l - j ) , (5.22a) g(m, c) = ac(l - c) - bmc(ccrit - c)e~ac . (5.22b) Similarly as in the previous case, I graph the dispersion function a, as seen in Figure 5.8. I note that I can dramatically increase the number of peaks n that form patterns by decreasing yu, i.e. the random motility of microglia. Figure 5.8: Dispersion relation for study Case 2 for several values of the bifurcation parameter a. Other parameters: a = 0.1, b = 1, c = 0.5, k = 1, Ccrit = 0.5, r = 5, x = 10, D=l . 3 Conclusion The extended Keller-Segel models have to undergo a more extensive investigation. I would like to apply all the mathematical tools used to study the Keller-Segel chemotactic system to these models since they show strong pattern formation. The form of the limiting kinetics function incorporates a couple of new features: microglia cells are able to either produce or remove the toxic chemical, depending on the model 82 0 2 4 6 8 10 12 q Figure 5.8: Dispersion relation for study Case 2 for several values of the bifurcation parameter fi. Other parameters:" a = 0.1, 6 = 1, c = 0.5, k = 1, Ccrit = 0.5, r = 5, x = 10, D = 1 -82a parameters. Also, I limited the chemical concentration increase, and implicitly, the mi-croglia density increase. My numerical simulations clearly indicate a 4-fold decrease in microglia density, compared to the Keller-Segel model that used a linear kinetics function. It is important to note that the Keller-Segel model and the extended Keller-Segel model that incorporates the limiting kinetics function have similar dispersion relations, i.e. q = 0 is the smallest value of the wavenumber for which the eigenvalue a is positive. However, the two systems develop different stationary distributions of microglia and chemical. Thus, numerical simulations generated for the extended Keller-Segel model with the limiting kinetics function indicate that the number of peaks formed in the chemotactic system initially correlates with the most unstable mode and later one or two large clusters of microglia and chemical form. Recall that the Keller-Segel model developed stationary peaks that correlated with the most unstable mode n. As yet, I do not know why the two models behave differently. I speculate this is due to the nonlinear terms that are neglected in the linear stability analysis, but play an important role in generating stationary distributions using numerical simulations (the limiting kinetics function is highly nonlinear). If I extend the Keller-Segel model by using multiple kinetics functions, then a different shape of the dispersion relation is obtained. In this case, the dispersion relation predicts the formation of many aggregates. This is due to the fact that the linear instability condition is satisfied for positive values of the wavenumber q, far away from zero. 83 Chapter 6 The Chemoattraction-Repulsion Model In this chapter, I introduce a more complex chemotactic model for the interactions be-tween one cell type and two chemicals, one attractant and one repellent. The motivations for studying this model are both mathematical and biological. Mathematically, the model leads to an interesting extension of traditional chemotaxis systems with new features, as proposed by A. Mogilner. Biologically, it is recognized that attraction and repulsion mechanisms are vital in the development of the central nervous system and neuronal growth cones (Mark et al. 1997; de Castro et al. 1999; Bagnard et al. 2000)). For example, different forms of netrins and semaphorins are either chemorepellents or chemoattractants for growth cones during development. In my biological system, I want to assess possible attraction and repulsion roles for glial cells. Thus, many soluble, diffusible substances (IL-1/5, IL-6, /3-amyloid and others) are secreted by and attracted to microglia and astrocyte cells (Mrak et al. 2000; Lee et al. 1993; Mackenzie et al. 1995; Stalder et al. 1999; Smits et al. 2000), but the mechanisms are not yet fully known from experiments. In addition, circumstantial evidence suggests a chemorepellent role for T N F - a (Chicoine et al. 1995). 84 1 Model Equations and Mathematical Tools I introduce a model of a single population of chemotactic cells and two types of chemicals that they secrete. This model is based on joint work with Alex Mogilner. In ID, the model equations are: dm d2m d ( dcx dc2\ - r n ^ ^ - d x [ X i m ^ ~ X 2 m ^ ) ' ( 6 - l a ) dc d2C' ~dt = Di~dx^ + a i m ~ b i C i ^ = 1 ' 2 > (6.1b) where Ci(x, t) are the chemical concentrations of attractant (i — 1) and repellent (i = 2), DiS are diffusion coefficients for the chemicals, \x is the cell random motility, X\iXi a r e chemotactic coefficients for attraction to chemical 1 and repulsion from chemical 2, and di, bis are rates of production and decay of the chemicals, respectively. For this model I will assume the same no flux boundary conditions: dm dx dcj dx = 0 , (6.2a) x=0,L = 0 . (6.2b) x-0,L This means that the cells and chemicals do not leave or enter the domain. The initial conditions are: m{x, 0) = u(x) , 0 < x < L , (6.3a) Ci(x, 0) = Wi(x) , 0 < x < L , (6.3b) where u(x) and Wi(x) are some sufficiently smooth functions. The mathematical tools that I use for analysis are as follows: • the method of Green's function solves the chemical equations (Zwillinger 1989). • linear stability analysis provides conditions for growth/decay of initial aggregates (Keller and Segel 1970; Edelstein-Keshet 1988). 85 • numerical simulations generate stationary non-homogeneous solutions. For this particular model, I used real, biological parameters in the simulations. 2 Nondimensionalization Investigation of the model given by Equations 6.1 starts with its nondimensionalization, in order to reduce the number of parameters (Murray 1993). I introduce the following notations: U = y / D J h , L 2 =^/D2~Jb~2, r =L22/p, (6.4) where L \ and L 2 represent, respectively, the effective ranges of the attractive and repulsive interactions between cells. Also, it will be seen that L 2 is a characteristic length scale of the problem and r is the time needed for a cell to move over this length scale. I then define the dimensionless quantities x t m x = , t* = - , m* = ^ , c* = ^ , (6.5) L 2 r m Ci where Ci = — m i = 1, 2 . (6.6) Because m and c\ are connected through relation 6.6, it follows that m = m , ci = ci , c 2 = c 2 , (6.7) where ( m , Ci,c" 2) represents the homogeneous steady state of the original System 6.1. Substitution of 6.5 into System 6.1 leads to the following nondimensional system of equations, after I drop the asterisks (see details in Appendix D l ) : dm d2m d dt dx2 dx . dci dc2 A l — - A 2 — m dx dx (6.8a) e2— = - ^ + r n - c 2 , (6.8c) 86 where , Xiairn . X2a2m p ii L2 Ai = — j — , A2=.—— , el = — , e2 = - - , a = -=-.. (6.9) / i O i fj,b2 JJi D2 Li The new boundary conditions are: dm dx dcj dx The new initial conditions are: x=0,-±-L2 = 0 , (6.10a) = 0 . (6.10b) m(x, 0) = u(x) , 0<x<^-, (6.11a) L2 Ci(x,0) = Wi(x) , 0<x<-^. (6.11b) L2 In order to continue my analysis, I consider the following simplification: I will assume about System 6.1 that chemicals diffuse more rapidly than cells, i.e. \i <C Di, which im-plies ti <C 1 in System 6.8. Following Cross and Hohenberg (1993), this condition implies that, on a fast time scale, chemicals diffuse very rapidly and reach a stationary distri-bution, while the cells do not undergo any movement. Therefore, a quasi-stationary distribution of chemicals develops. Such a simplification will allow me to reduce the system to the "shadow system", a single equation for m, by finding Green's function solutions for C\ and c2. I should also note that, using the scales in 6.5, the homogeneous steady state of the nondimensional System 6.8 is (in, c\, c2) = (1,1,1). 3 The Method of Green's Function By assuming a quasi-stationary distribution of chemicals, I obtain the following equations for the chemicals Q d2Ci 0 = -Q-J - OJjCj + axm , i = l,2, (6.12) 87 where ct\ = a 2 and a 2 = 1- Equation 6.12 is a second order, linear, non-homogeneous ordinary differential equation in Cj. I recall from Chapter 3, Section 6.1 that it can be solved using the method of Green's f u n c t i o n on an infinite domain (—00, 00) (Zwillinger 1989). Since I am concerned with a finite domain [0, L], I look for an additional assump-tion that makes possible the use of Green's function on an infinite domain. The condition is related to the effect of the boundaries on the chemical distribution. Let us analyze the following dimensionless ratio (Lee et al. 2001): where D i and b{ are the diffusion and decay rates of chemicals, and L is the domain length. Equation 6.13 can be interpreted as the ratio of the chemical half-life time to the char-acteristic time for the chemical to diffuse over the domain length. In Chapter 7, I show that the following sets of values are reasonable: b\ = 0.18 / m i n , D \ = 900 pm2/min, 62 = 0.045 / m i n , D 2 = 900 pm2/min, and L = 700 jim. These parameter values lead to ri = 0.01 and r2 = 0.04. It follows that r ; < 1. The interpretation of this last in-equality is that the chemical decays faster than it diffuses across the domain. Therefore, the boundary conditions have little influence on the dynamics of the c h e m i c a l (Lee et al. Having found this important condition regarding the effect of the boundary conditions on the chemical distribution, I can now reason that the desired finite domain Green's function is just a small perturbation of the function defined for the infinite domain (Morse and Feshbach 1953). The details of the chemical distribution derivation are shown in Appendix D2. The desired solution can be expressed as where fi(x) = aifn(x) and Gi(x,x') is the associated Green's function, given by (6.13) 2001). (6.14) Gi(x,x') 1 e-%/a~i\x-x'\ (6.15) 88 and shown in Figure 6.1. The Green's function represents the chemical concentration at 0.4 -3 -2 Figure 6.1: Green's function G(x,x'). x resulting from the production of this chemical by cells at point x'. The solution, q , of Equation 6.12 reads The form of the chemical Equation 6.16 indicates that at a point x, the chemical con-centration Ci(x) is found by superimposing concentrations resulting from cells secreting at all points x', —oo < x' < oo. Therefore, I can say that the interactions between the cells are non-local. I will use Equation 6.16 to reduce the complexity of the model to a single equation for m(x, t). This is done in the following section. I should mention here that, if I were to solve the dimensional chemical Equation 6.1b, then the decaying Green's function solution, with source at x' = 0, would have been given by Go/e, i.e. for x = Li, Gi(x,0) has fallen to 1/e of its value at x = 0. Thus, I choose (6.16) (6.17) where Go is the value of Gi(x, 0) at the origin. It follows that for x = y/Di/bi, Gi(x, 0) = 89 x = L 2 to be the characteristic length scale of the problem. This fact is visualized in Figure 6.2. Figure 6.2: Plot of the dimensional decaying Green's function G(x,0). For x = L2, the Green's function decays to 1/e of its value at the origin. L2 represents the characteristic length scale of the model. 4 The Shadow System I turn my attention to the cell density Equation 6.8a. Using Equation 6.16, the cell equation can be expressed, as shown in Appendix D3, as an integro-partial differential equation of the form dm d2m d dt dx2 dx (vm) (6.18) where A 2 A 2 f ° ° v — —(K * m) = — / K(x — x')m(x') dx' , 2 2 J - o o (6.19) represents non-local drift velocity, and the function K(x) is the convolution kernel given by K(x) = sign(—x) (Ae -o\A _ P-V (6.20) 90 The kernel K(x) is an odd function and it represents the effect of a source at x = 0 on the motion of a cell at the point x. Here the parameters have been absorbed into X i ai D2 L2 , v . A = — , a = — . (6.21) X2 0-2 Dx Lx A and a are dimensionless parameters. A represents the ratio of effective strengths of the attraction and repulsion, while a represents the ratio of the spatial ranges of the repulsion and attraction. The model given by Equation 6.18 is a spatially non-local integro-differential model. This means that each cell interacts with many other cells. These interactions could be of two types: attraction or repulsion, and they can be visualized when plotting the kernel 6.20 (see Figure 6.7). This important feature of the model will be also analyzed in the parameter space Aa shown in Figure 6.5. Equation 6.18 is called the shadow system since it represents a simplification of the full System 6.8. In the following section, it will be seen that linear stability analyses of both systems, i.e. the full System 6.8 and the shadow System 6.18, lead to the same results. 5 Linear Stability Analysis 5.1 The Full System I consider spatially non-homogeneous small perturbations m! and c- of the homogeneous steady state (1,1,1) of the full System 6.8 such that m(x,t) =1 + rri = 1 + Meatcos(qx) ,• ' (6.22a) Ci{x, t) = 1 + cj = 1 + Cieat cos(qx) , (6.22b) where • q = mr/L is the wavenumber and the positive integer n is the mode. 91 • a, the eigenvalue, represents the growth rate of the perturbations. Its sign will determine whether I'll obtain growth or decay of perturbations. • M and Ci are amplitudes of the perturbations at time t = 0. I linearize the full system by substituting Equations 6.22 into System 6.8 and retaining only the linear terms. I obtain the following linear system (after dropping the primes): (6.23a) e2 dm ~dt dci ~dt dc2 ~dt d2m dx2 d2cx dx2 d2c2 dx2 dx2 dx2 + a2(m — c i ) , m - c2 . (6.23b) (6.23c) Further substitution of the perturbations shown in Equations 6.22 leads to a system of three equations with the unknowns M, C\ and C2 that can be written in matrix form as: fa \ V f M \ Ci = 0 . (6.24) 2 - A x q 2 A 2 q 2 a2 - (a 2 + q2 + exa) 0 1 0 - (1 + g 2 + e2cr) J System 6.24 has a trivial solution ( M , Ci,C2) = (0,0,0). A non-trivial solution exists if and only if the determinant of the coefficients matrix is zero, i.e. fa - A x q 2 A 2 q 2 det a 2 - (a 2 + q2 + exa) 0 V 1 0 - ( l + g 2 + e2cr) ) Calculation of the determinant leads to a cubic equation in a: a3 + aa2 + Pa + 7 = 0 , where = 0 (6.25) a = Q2 + P = " 1 . -1 7 = t i t 2 ~2 1 „2\ 1 £2 (6.26) (6.27a) 1 62 A2 _ a 2 A x e2 ei eie2 ; i + g 2)(a 2 + g2) ,(6.27b) i)} Q2 . (6.27c) 92 Equation 6.26 represents the dispersion relation associated with the full chemoattraction-repulsion System 6.8. I know that instability of the system occurs if growth of perturba-tions develops, i.e. Real(a) > 0. To analyze the signs of the cubic roots of Equation 6.26, I make the following simple observations: • if O i , (T2 and <73 are the three roots of Equation 6.26, then a = - ( o i + o2 + , (6.28a) 8 = <7i<72 + CT1CT3 + 0203 , (6.28b) 7 = - ( C T I C T 2 C T 3 ) . (6.28c) • if o i , o"2 and cr3 are all negative, or o\ < 0 and 0-2,3 £ C, then 7 > 0. • if one root becomes positive, then 7 < 0. Therefore, the linear instability condition for the full System 6.8 is 7 < 0 . (6.29) Substitution of Equation 6.27c into the instability condition 6.29 gives the aggregation condition — — > — . (6.30) I will immediately see that the same aggregation condition is obtained when I do a linear stability analysis for the reduced system, i.e. the shadow system 6.18. 5.2 The Shadow System Similarly, I consider small perturbations of the homogeneous steady state rh = 1, i.e. m(x, t) = 1+m', where m! = M e a t cos(qx). Instability develops if growth of perturbations occurs, i.e. Real(cr) > 0. Substitution of the perturbation relation into the shadow 93 System 6.18 gives dm' dm' d dt dx2 dx A 2 (K * (1 + m'))(l + rri) (6.31) I want to linearize Equation 6.31, and therefore I only retain linear terms that appear in this relation. Mathematical calculations, as shown in Appendix D4, lead to dm' d2m' d dt dx2 dx Furthermore, substitution of the perturbation ml = M e a t cos(gx) into 6.32 gives ( A 1 o M e a t cos(qx) = - q 2 M e a t cos(gx) + A 2 q 2 M e a t cos(gx) (6.32) (6.33) g 2 + a2 g 2 + 1 Division of Equation 6.33 by M e a t cos(gx) yields the linear dispersion relation for the shadow system: A 1 N (6.34) a = - g + g A 2 \qz + a1 qz + 1 The aggregation condition associated with the shadow System 6.18 is then A 1 1 > (6.35) g"z + a z g z + 1 " A 2 Thus, when comparing the two aggregation conditions 6.30 and 6.35, I observe that they are identical. I can therefore continue my study of pattern formation of the full chemoattraction-repulsion System 6.8 using the reduced shadow System 6.18. 5.3 T h e Bi fu rca t ion Func t ion I introduce the bifurcation function A = 2 , 2 ql + a2 g 2 + l ' (6.36) so that, from Equation 6.34, a(g)=g^ (A 2 a*(g) -1) (6.37) 94 For a better understanding of the important role the bifurcation function 6.36 plays in the linear stability analysis, I carry on with a detailed description of cr*(q). Inspection of the Relation 6.36 reveals that a* is expressed as a difference of two ratios. Let (6.38) (6.39) q c. / i and fa for A = 2, a = 1.5. Figure 6.3: Graphs of the functions fa and fa for different values of A and a. (a) a* < 0; (b) a* > 0 and has a maximum at q = 0; (c) a* has a maximum at q ^ 0. For different values of A and a, these two ratios can be plotted together to illustrate the sign of a* and its properties. Thus, in Figure 6.3a, a* is always negative, and therefore 95 no pattern formation can occur. Figure 6.3b shows that a* is always positive and has a local maximum at q = 0. Finally, in Figure 6.3c I see that a* changes sign. Let qr be the zero of the function a*(q), i.e. <J*(qr) = 0. Then l a 2 - A . A , The existence of real valued qr depends on A and a, and requires that (a 2 — A ) / ( A — 1). 1 1 1 1 1 1/^ =1, "max / / .J. /q min / Qmax - / / / -3 1 • 2 3 4 5 € Figure 6.4: The bifurcation function a* for A = 2, a = 1.5, and the line l /A? for different parameter sets. For a* < I/A2, the homogeneous steady state is stable. For a* = 1/A2, there is bifurcation to instability. For a* > I/A2, aggregation occurs. Recall from Inequality 6.35 that whenever a*(q) > l / A 2 , pattern formation will develop. It follows that additional information about the parameters A and a can be obtained by analyzing how the curve a* and the line l / A 2 relate to each other graphically. Thus, I plot the two curves for different values of the parameters A and a, as shown in Figure 6.4. In this figure, we observe that for 1/A2 = h, a* < 1/A2, and the aggregation condition 6.35 is not satisfied. It follows that a is negative and therefore the homogeneous steady state remains stable. The line I/A2 = l2 is tangent to the curve a* and their intersection represents the bifurcation to growth of perturbations. The bifurcation occurs for a critical 96 wavenumber qc for which o~*(qc) = cr^ a x. The critical wavenumber qc is given by 'a2 - VA V A - I ' whenever the quantity under the square root is positive. Then cr^ a x is (6.41) (6.42) As the line 1/A 2 moves downward, more values of the wavenumber q can be found for which cr*(q) > l / A 2 . The aggregation inequality 6.35 is then satisfied and pattern formation should occur. In this case, if 1/A 2 = h for example, there is an interval of wavenumbers (qmin, Qmax) for which aggregates will form. 5.4 T h e Parameter Space Aa In the previous section we have seen that it is convenient to analyze the sign of a* in the parameter space A a . The goal is to explore what values of a and A give pattern formation. Therefore, I plot four important curves in the parameter space A a , namely A = 1, a = 1, a = A 1 ! 2 , and a = A 1 / 4 , as shown in Figure 6.5. Region 1 No pattern formation Region 3 Periodic pattern formation Figure 6.5: The parameter space Aa. In regions 2 and 3, whenever a* > 1/A 2, aggregation of cells occurs. In region 1, situated to the left of the curve A = 1 and above a = A 1 ! 2 , the bifurcation function a* is negative for all values of the wavenumber, as seen in Figure 6.6a. The 97 'X, \ \ \ \ Oc = ° a. Region 1: A = 0.5, a = 0.9. b. Region 2: A = 2, a = 0.5. c. Region 3: A = 2, a — 1.5. Figure 6.6: The bifurcation function a* in different regions of the parameter space Aa. In case a, a* < 0 implies that the uniform state is stable. In b, there is a tendency to form large (low wavenumber) clusters. In c, the fact that a* has a local maximum at qc » 2 implies that patterns with this wavenumber will be first to form when bifurcation occurs. 98 X a. Region la: A = 0.5, a = 1.5. x b. Region lb: A = 0.5, a = 0.8. e. Region 2e: A = 20, a = 2. f. Region 3: A = 1.5, a = 1.75. Figure 6.7: The kernel K(x) in different regions of the parameter space Aa. This function represents the effect of a cell (which secretes chemicals) at x = 0 on the velocity of other cells in its neighborhood. K > 0 rep-resents motion to the right, K < 0, motion to the left. (See also the corresponding arrows along the x axis). 99 homogeneous steady state is stable and therefore no pattern can form. In region la , i.e. above the line a = 1, the interaction between cells is of pure repulsion (Figure 6.7a). In region lb , when All2 < a < 1, there is local repulsion between cells and weak attraction at larger distances (Figure 6.7b). In region 2, situated below the curves a = All2 and a = Al'\ when a* > l/A2 for different values of q, instability occurs (see Figure 6.6b). In this region, the loss of stability occurs at qc = 0 and a* has an absolute positive maximum at this critical point. I therefore expect large scale pattern formation to develop (Cross and Hohenberg 1993). In region 2c, when A < 1, there is local repulsion between cells and strong attraction at larger distances (Figure 6.7c). I can also have pure attraction between cells, as in region 2d (Figure 6.7d). Finally, in region 2e, when a > 1, there is attraction at small distances and weak repulsion at larger distances (Figure 6.7e). In region 3, situated to the right of A = 1 and above a = Al/A, a* > 1/A2 for larger values of q (see Figure 6.6c). The bifurcation function has an absolute positive maximum at qc. Loss of stability occurs at qc ^ 0, and therefore periodic pattern is expected to develop. In this region there is only local attraction between the cells and repulsion at larger distances, as seen in Figure 6.7f. 6 Biological Interpretations of the Aggregation Con-dition Similar to my discussion regarding the Keller-Segel aggregation condition, I now examine and interpret biologically the aggregation condition given by Inequality 6.35. There are many parameters in the linear instability condition, and therefore it is more difficult to biologically interpret the inequality. To possibly suggest targets for therapy intervention in Alzheimer's disease, I should look for conditions such that no plaques 100 would form in the tissue. Thus, from Figure 6.5, I observe that region 1 is of interest to us since there is no pattern formation there. In particular, if the following two inequalites are simultaneously satisfied, i.e. a > 1 and A < 1, then clearly no plaques will form. These two inequalities make the left-hand side of the aggregation Inequality 6.35 negative. To interpret biologically what a > 1 and A < 1 mean, I have to analyze the definitions of both parameters. I have: repulsion spatial range . t „ . a = > 1 , (6.43a) attraction spatial range attraction strength . A = — < 1 . 6.43b) repulsion strength Mathematically, this is equivalent to: L 2 > L i , (6.44a) X i a i L - l b i < X2a2L\bx . (6.44b) These two inequalities suggest that no plaques form in the tissue if the following condi-tions are satisfied by the system: • the repulsion spatial range should be L A R G E R than the attraction spatial range. • since L{ = y / D i / b i , I propose that the decay rate 62 of the chemorepellent should be smaller than the decay rate b\ of the chemoattractant. • increase the production rate a 2 of the chemorepellent. Thus, one should thoroughly investigate the role of chemorepellents in the brain. So far, only TNF-ct has been identified as a possible chemorepellent under certain conditions (Chicoine et al. 1995). I suggest that a decrease of the decay rate and an increase of the production rate of chemorepellents in the brain could be beneficial in controlling plaque formation in Alzheimer's disease. 101 7 Numerical Simulations The numerical simulations shown in this section were generated for the nondimensional chemotactic model given by Equations 6.8, the boundary conditions 6.10 and random initial conditions, using real biological parameters. I used the same numerical scheme as described in Chapter 4, Section 7. The derivation of the parameters requires a separate, detailed discussion and constitutes the subject of the next chapter. I used two sets of parameters, and their values are given in Table 6.1. Using these parameters, I graph the bifurcation function cr*, as seen in Figure 6.8, in order to find the range of the wavenumber q where pattern formation could occur. For the parameter set 1, I find that for 5 < n < 9, aggregation of cells should develop. For the parameter set 2, aggregation should occur for 4 < n < 21. For both sets of parameters, the most unstable wavenumber corresponds to n = 6. Parameter Set 1 Set 2 Al 37.50 150 A2 100 400 ei 0.0367 0.0367 0.0367 0.0367 a 2 2 A 1.5 1.5 domain length 5 5 Table 6.1: Dimensionless parameter values for the chemoattraction-repulsion .model Equations 6.8. I discuss now the difference between the numerical results obtained using the two sets of parameters. In Figures 6.9f and 6.10f, I plotted the stationary distributions of microglia and chemicals, and I observe that 3 aggregates, and respectively 4.5 aggregates, formed 102 0 005 0.0025 0 0 2 4 6 8 10 12 14 q Figure 6.8: The bifurcation function a* for A = 1.5, a = 2, and the line y = 1/A2 for A 2 = 100 and A 2 — 400. For a* > l/A2, aggregation occurs. in the domain. From Table 6.1, I see that the nondimensional parameter A2 has different values for the two sets of parameters. This leads to two different values of the line y = 1/A2, as shown in Figure 6.8. When A2 = 100, the values of q G (qmm, ?max) on the line y = l/A2 are closer to the bifurcation wavenumber qc, and therefore it is expected that the system will become unstable and form clusters with wavelength given by qc- When A2 = 400, the line y = l/A2 drops and qmax increases in value. From Segel (1984), it is known that in such a case, the number of peaks forming in the domain is approximately correlated with g m a x . My numerical simulations clearly reflect these last considerations: in Figure 6.9f there are 3 peaks (most unstable mode n = 6); in Figure 6.10f there are 4.5 peaks ( n m a x = 21). An additional note has to be made regarding the number of peaks associated with the wavenumber g m a x > qc. The number of peaks that I observe in Figure 6.10f is due to effects of nonlinearity (Segel 1984). For A2 = 400, 4 < n < 21, and therefore the system can evolve to any mode n in this interval. If I approximate the random initial conditions as resulting from superimposing several periodic modes for the different values of n G [4, 21], then it is difficult to indicate what mode will grow into a stationary pattern. Clearly, the i 1 1 1 1 r 103 (a) t=0. (b) t=0.02. microglia / \ attractant j I \ repelent \ y (c) t=0.04. (d) t=0.14. E 005 |l rnjcraglia ,y i \ ! 1 a«ractant 1 1 1 1 j \ | | / repellenl | \ 1 1 i I / \ i i / \ J \ l\ \j \ i \ 1 \ V / J \ (e) t=l . (f) t=2. Figure 6.9: Spatial distributions of microglia and chemicals for the nondimensional chemotactic model, using random initial conditions. Parameters: A\ = 37.5,^2 = 100, = 0.0367, a = 2 , L / L 2 = 5,e = 0.002. Homogeneous steady state: m = Ci = 1. Random initial conditions: m(x,0) = rh + ea;ran, Ci(x, 0) = Ci + ex r a n . 104 • A / 1 1 1 ^ / 1 attractant/ \ 1 ' •\ | 1 X \ I - \ _jL-4. • —' 1 '1 microglia 1 | repellent 1 1 • { \l V . v v 1/ 'j 51/ (b) t=0.01. j 1 microglia/ I 1/ / i r / \ / \ / ' \ repellent \' \ A I M / (c) t=0.03. ,1 microglia « II— II A j I repellent j \ j jattractant i \ I \ J (d) t=0.19. .^microglia repellent E 00025 1L-microglia j j | repellent 1 I \ attractant 1 ir 11 v \ ; i A V \ (e) t=0.2. (f) t=l . Figure 6.10: Spatial distributions of microglia and chemicals for the nondimen-sional chemotactic model, using random initial conditions. Parame-ters: Ai = 150, A2 = 400, et = 0.0367, a = 2,L/L2 = 5,e = 0.002. Homogeneous steady state: TO = Ci = 1. Random initial conditions: m(x, 0) = TO + ea;r a n, Ci(x, 0) = Ci + erc r a n. 105 nonlinear terms that I neglected in my linear stability analysis become important, and therefore a n o n l i n e a r analysis should be developed. 8 Conclusion In conclusion, analysis of the chemoattraction-repulsion model predicts that periodic pattern formation occurs for local attraction and long range repulsion interactions be-tween cells and chemicals. I should emphasize that the linear stability analysis is local to the homogeneous steady state and therefore I cannot anticipate the behaviour of the chemotactic system far from the homogeneous steady state. Also, the periodic pat-tern formation predicted in region 3 of the parameter space A a occurs for values of the wavenumber in the interval (qmin, qma,x), as shown in Figure 6.4. Biological interpretations of the aggregation condition showed that an increase in the production rate of the repellent might prove to be beneficial in slowing the progression of plaque formation in Alzheimer's disease. Numerical simulations show that the number of clusters containing microglia cells and two types of chemicals correlates with the most unstable wavenumber whenever q G ( t 7 m i n , qma.x) is close to the bifurcation value qc. Nonlinear effects influence the dynamics of the system when the wavenumber i 2 m a x becomes much bigger than the bifurcation wavenumber qc. An additional accomplishment of the numerical simulations is the use of parameter values that are biologically relevant to Alzheimer's disease. In the next chapter, I will show how I derived the parameters, and I will make the connection between the theoretical model and the biology of Alzheimer's disease. 106 Chapter 7 Applicability of Models to Alzheimer's Disease A final goal of the thesis is to include biological data, available in the literature, into the models. Specifically, I would like to identify and estimate parameters that appear in the models; incorporate them into the models to generate numerical simulations; compare theoretical predictions with real biology and, consequently, indicate and discuss the limitations of the models. The parameters that I was able to gather are shown in tables in the first section. They were divided into three main categories: parameters that are currently available in the literature and therefore have known values; parameters that can be rigorously calculated; finally, parameters that are estimated. I display the parameters in groups based on the function they play in the models, e.g. diffusion coefficients, production rates, etc. Included are certain related parameters, not used directly in the models, that might prove useful in future developments of the models. 1 Parameter Derivation In Table 7.1, I show the molecular weight of proteins and cytokines. These constants will be used to estimate other parameters that are known to have particular values based on their specific molecular weight. For example, the effective diffusivity of molecules 107 can be estimated using molecular weight. These diffusion coefficients are given in Table 7.2. Molecule Type Value Source ApoE 37 k D a Krul and Tang (1992) A P P 125 k D a Dickson et al. (1995) /3-amyloid 4 k D a Jordan et al. (1998) IL-1/? 17 k D a Benveniste (1995) IL-6 26 k D a Benveniste (1995) T N F - a 17 k D a Benveniste (1995) Table 7.1: Molecular weight of proteins and cytokines. Molecule Type Value Comments Reference 6-amyloid 1500 ^ mm C a l c u l a t e d Goodhill (1997) IL-1/3 900 mm C a l c u l a t e d Goodhill (1997) IL-6 810 mm C a l c u l a t e d Goodhill (1997) Moghe et al. (1995) T N F - a 900 mm C a l c u l a t e d Goodhill (1997) Table 7.2: Effective diffusivity coefficients of various molecules. I calculated the effective diffusivity of molecules as follows: 1. The diffusion coefficient of a molecule in cytoplasm scales approximately as the inverse of the cubic root of its molecular weight (Goodhill 1997). 2. Experiments show that a molecule with a molecular weight of 0.5 kDa has a diffusion coefficient of approximately D = 10~6 cm2/s (Goodhill 1997). This result helps us to scale any diffusion coefficient. I call this particular molecule the scaling molecule. 108 3. I use the fact that 10 6 cm2/s = 6000 pm2/min . 4. Thus, I can calculate the diffusion coefficient of a molecule if I know its molecular weight. For example, IL-1/? has a molecular weight of 17 kDa and its diffusion coefficient is given by: Mwt i 17 = — = 34 Mwt 0.5 1 _ 1 1 ^ M w t i ~ y/U v 'Mwt £>i = 0.3 D Dx = 1800 ^ , mm where Mwt i , D\ and Mwt, D are the molecular weights and diffusion coefficients of IL-1/3 and the scaling molecule, respectively. 5. The calculated diffusion coefficient is not accurate because my molecules are not diffusing in the cytoplasm. I have to incorporate a correction factor of one half in the above calculations due to the tortuosity of the neuronal extracellular space (Sykova 1997; Mazel et al. 1998; Sykova et al. 1998; Nicholson and Sykova 1998). I then obtain the effective diffusion coefficients shown in Table 7.2. In Table 7.3, kinetics parameters, such as size, speed and concentrations in various mediums are given for different cell types. Most of the values were found in the literature, except for the microglia speed which was estimated using the macrophage speed. Tables 7.4 and 7.5 present several parameters for cytokine receptors. Most of the values were found in the literature and will be used for future calculations. The receptor concentration parameter rc, shown in Table 7.4, was calculated as follows: 1. Let V be the volume of a cell with membrane receptors. Here I am interested in the volume of a microglia. 109 Cell Type Description Value Comments Source Microglia Size (cell body) 10-15 jim Human (H) Streit (1995a) Microglia Size (reactive) 90-100 um Human Streit (1995a) Astrocytes Size (reactive) 15-25 pm Human Marzolo et al. (2000) Macrophage Speed 2 fim min Farrell et al. (1990) Microglia Speed min Estimated N / A Astrocytes Speed 0.1-0.2 ^ mm Kornyei et al. (2000) Microglia Density 100-350 cells W4fj,m2 H. hippoc. (AD) Itagaki et al. (1989) Microglia Density o.5 cells l O V m 2 H. temporal neocor. (AD) Mackenzie et al. (1995) Table 7.3: Cell kinetics parameters for different cell types. 2. A normal microglia is roughly 10x10x10 \im. This gives an estimate for the volume of a microglia to be V « 1000 pm3. 3. Let r be the number of receptors per cell as given in the second column of Table 7.4. 4. Since I would like to find the receptor concentration in units of n M (nano mole), I have to express the number of receptors per volume in this unit. The conversion Cytokine Receptors (r) / Cell Receptor Concentration (r c) Source IL-1/3 5000/macrophage 8.3 n M Benjamin et al. (1990) IL-6 500/macrophage 0.83 n M Yamaguchi et al. (1992) T N F - a 4000/macrophage 6.64 n M Michishita et al. (1990) Table 7.4: Receptor kinetics parameters for cytokines. 110 Cytokine [k-]= 1 min [ * + ] = I nM—min nM Source IL-1/3 0.004-0.04 0.004-0.04 1 Benjamin et al. (1990) IL-6 0.047 0.012 4 Hammacher et al. (1996) T N F - a 0.01 0.0036-0.05 0.2-2.8 Michishita et al. (1990) Pennica et al. (1992) Table 7.5: Binding rates (k+), unbinding rates (&;_) and dissociation rates [KD) for cytokine receptors. is done using Avogadro's number. One mole consists of 6.02 x 10 2 3 molecules Cytokine Produced by Value Comments IL-1/3 Microglia (stimulated by LPS) 300 ^ 0 . 2 x l 06 cells ml for Ah IL-6 Microglia (stimulated by LPS) 2000 ^ 0 . 2 x l 06 cells ml for Ah T N F - a Microglia (stimulated by LPS) 400 & ml 0 . 2 x l 0 6 cells ml for Ah IL-6 Astrocytes (stimulated by IL-1/3) 300 & ml 0 . 2 x l 0 6 cells ml for Sh T N F - a Astrocytes (stimulated by IL-1/?) 80 & ml 0 . 2 x l 0 6 cells ml for 8h Table 7.6: Concentration of cytokines produced by microglia or astrocytes, in the presence of a stimulating factor. Al l values from Lee et al. (1993). I l l (receptors in this case). I obtain: 1 r 1 i 1 r , i rc = r receptors x — x —— — \mole\ H V[fim3i 6.02 x 10 2 3 receptors1 J r r 1 , 1 rw 7-n ^ T i r M 7 7 * n - , ^ Q r n - ^ \ = v f e 1 x 6mVW*[M x L] x 1 0 [ — ] x 1 0 V 1 = ^ x 1.66 nM . 5. In the previous conversion I used: • 1 mole = 1 M x 1 L. • 1 L = 1 d m 3 = 10 1 5 / u m 3 . • 1 M = 109 n M . 6. For example, I find the receptor concentration of IL-1/3 on a microglia cell to be = x i.66 nM = 8.3 nM . (7.1) 1000 v ' Table 7.6 contains only known parameters. They represent the concentration of cy-tokines given in pg/ml (picogram per milliliter) produced by microglia or astrocytes, in the presence of a stimulating factor. Using these concentration values, I can derive production rates of cytokines. I deduce the production rates in Table 7.7, using the following calculations: 1. If there are n cells secreting for a period of x hours, then the amount of chemical produced by one cell per minute will be: I F T X — x - - = — — — — - M - ? - per cell . ml n cells x h 60 x n x x mm 2. For example; I find that one microglia cell secretes the following amount of IL-1/3 per minute: 3 0 0 r P9 , _ 6 . 2 5 x 1 0 - e r P9 , 60 x 0.2 x 106 x 4 L m m J ' L m m J 112 Cytokine Produced by One Cell Value IL-1/? Microglia (stimulated by LPS) 6.25 x l O " 6 M-mm IL-6 Microglia (stimulated by LPS) 41.67 x l O - 6 -22-mm T N F - a Microglia (stimulated by LPS) 8.33 x ICT 6 -E2-mm IL-6 Astrocyte (stimulated by IL-1/?) 3.125 x ICT 6 mm T N F - a Astrocyte (stimulated by IL-1/?) 0.833 x 10" 6 mm Table 7.7: Production rates of cytokines produced by microglia or astrocytes, in the presence of a stimulating factor. Al l values were calculated from Lee et al. (1993). In Table 7.8, all cytokine decay rates were calculated. In general, chemicals decay because cells within the brain absorb them through receptors. Since the available liter-ature does not supply this type of information, I made the following assumptions and calculations: Cytokine Absorbed by Value IL-1/? Macrophage 0.033-0.33 -K-mm IL-6 Macrophage 0.01 -X-mm T N F - a Macrophage 0.024-0.33 -X-mm Table 7.8: Decay rates of cytokines. 113 Parameter Descr ip t ion C e l l T ype Value Source X Chemotaxis towards IL-8 Human neutrophils 6 780 ^ m 2 Moghe et al. (1995) Random motility Human neutrophils 33 ^ mm Moghe et al. (1995) Table 7.9: Additional parameters for cell movement. 1. The absorbing cell has receptors on its surface. I previously calculated the receptor concentration, r c , in units of nM (Table 7.4). 2. A chemical binds to an absorbing cell with a binding rate k+ (Table 7.5). 3. I therefore define the decay rate as b = k+ x rc [1/min]. In Table 7.9, I give two additional parameter values for microglia which were directly found in the literature. However, the random motility coefficient fj, was corrected to account for the tortuosity of the neuronal extracellular space. 2 Applicability of the Chemoattraction-Repulsion Model Throughout the thesis, I proposed several mathematical models in order to theoretically investigate the dynamics of chemotactic systems. Some of the analytical results were verified by numerical simulations. It is important to note that all parameter values used in these simulations were arbitrarily chosen to satisfy the conditions and assumptions derived from the models, except in the case of the chemoattraction-repulsion model where real, biological parameters were used. In the previous section, I have shown how I derived all the parameters; I can now investigate, for the chemoattraction-repulsion model, whether the parameters found in the literature would be consistent with the 114 development of aggregates, e.g. as seen in Alzheimer's disease. 2.1 Model Parameters The chemoattraction-repulsion model was developed to study the interactions between one cell type and two chemicals. More specifically, one chemical is a chemoattractant for the cell, while the other chemical is a chemorepellent for the cell. Microglia are the cells represented in this model. I speculate that the chemoattractant chemical is IL-1/3 and the chemorepellent chemical is T N F - a . As already mentioned, the IL-1/5 cytokine plays several roles in the initiation of plaque formation (Griffin et al. 1998): it activates microglia, inducing chemotaxis; it promotes production of T N F - a and IL-1/3 by microglia. Chicoine et al. (1995) found T N F - a to play a chemorepellent role in the brain. In Table 7.10, I present the dimensional parameter values that appear in the chemo-attraction-repulsion Equations 6.1. My numerical simulations were generated for the nondimensional system given by Equations 6.8. The parameters are: A i — —7—, A2 — —r—•) {<•*) ph 6 L = A ' € 2 = ^ ' ( ? - 3 ) X i a i D 2 L 2 , . A = -=-, a=j-, (7-4) X2 0,2 Di L x Li = VDi/h , L2 = ^D2/b2 , (7.5) where m is the dimensional microglia homogeneous steady state of the original System 6.1, and L i and L 2 represent the effective ranges of the attractive and repulsive interactions between cells. A represents the ratio of effective strengths (amplitudes) of the attraction and repulsion, while a represents the ratio of the spatial ranges of the repulsion and attraction. A l l the other parameters are given in Table 7.10. 115 Parameter Descr ip t ion Value Comments Microglia random motility 33 ^ min Known value from Moghe et al. (1995) Xi Chemoattraction 6 780 mm—nM Known value from Moghe et al. (1995) X2 Chemorepulsion ? Not Available D l IL-1/3 diffusion 900 mm Calculated using Goodhill (1997) D 2 T N F - a diffusion 900 mm Calculated using Goodhill (1997) ai IL-1/3 production rate (per microglia cell) 6.25 x 10~6 -22. mm Calculated using Lee et al. (1993) 0-2 T N F - a production rate (per microglia cell) 8.33 x 10~6 -22-mm Calculated using Lee et al. (1993) h IL-1/5 decay rate 0.033 - 0.33 -K-mm Calculated using b = k+ x rc b2 T N F - a decay rate 0.024-0.33 -X-mm Calculated using b = k+ x rc Table 7.10: Parameter values for the chemoattraction-repulsion model Equa-tions 6.1. The first step in my analysis is to verify whether the biological parameters satisfy the conditions used in the model. One assumption was that chemicals should diffuse more rapidly than the cell, i.e. e, <C 1. Using the values from Table 7.10, I calculate e'i « 0.0367 < 1 , e2 ~ 0.0367 < 1 . 116 Thus, this condition is satisfied for the biological parameters. A second condition was utilized when deriving the chemicals equations with the method of Green's function. Using the parameter values in Table 7.10, I have shown in Chapter 6, Section 3, that the ratios in Equation 6.13 are much less than one, and therefore the imposed assumption is satisfied. The next step in my investigations is to determine all the nondimensional parameters: Ai, A 2 , A , a. I find a using Equation 7.4 and values associated with IL-1/3 and T N F - a (see Table 7.10). Thus, a can be less than or bigger than 1, depending on the parameter values chosen for b± and b2, which implies that any of the three situation can arise: no pattern formation, large scale pattern formation and periodic pattern formation (from Figure 6.5). I focus my attention on the latter case, investigating the range of values that would lead to a > 1, i.e the repulsion range is greater than the attraction range between cells. Since D \ — D 2 , it follows that I can get a > 1 if b\ > b2, i.e. if IL-1/3 has a higher decay rate than T N F - a . Both decay rates have a range of values, and therefore I must decide on one value in order to make any progress. From Benjamin et al. (1990), I will assume that it is best to choose an average decay rate for IL-1/3, i.e. bi = 0.18/mm. Analysis of data available for T N F - a (Ding et al. 1989; Michishita et al. 1990; Pennica et al. 1992) suggests that an estimate for the decay rate for T N F - a is b2 = 0.045/'min. I then obtain a = 2. Unfortunately, I cannot calculate A directly since the chemorepulsion parameter X2 is not available. In region 3 of the parameter space A a , I know that 1 < A < a 4 (see Figure 6.5). Using this inequality, I decide to choose A = 1.5. In addition, I can predict that, if 1.33 < — < 21.33 , (7.6) X2 117 then periodic peaks should form. Inequality 7.6 is equivalent to 0-75 xi > X2 > 0.047 xi • (7-7) I know that Xi = 6 — 780 iirr? jmin — n M . It follows that 4.5 - 585 > X 2 > 0.282 - 36.66 . (7.8) I can therefore state that if the relative chemotactic coefficients associated with microglia and gradients of IL-1/3 (xi) and T N F - a (%2) are in the range 0.047 Xi < X2 < 0.75 Xi, then spontaneous instability leading to periodic pattern should develop. A = 1.5 gives the value of the ratio X1/X2 = 2, and choosing Xi = 6 pm2/min — n M , I obtain X2 = 3 11m2/min — nM. q0 = 3.5 0 2 4 6 8 10 12 14 q Figure 7.1: The bifurcation function a* for A = 1.5, a = 2, and the line 1/A2 for A i = 100 and A 2 = 400. For a* > l / A 2 , aggregation occurs. Analysis developed in the previous chapter showed that the linear instability condition for the chemoattraction-repulsion model is a* > l / A 2 , where a* is the bifurcation function given in Equation 6.36. Therefore, I search for values of A 2 that satisfy this inequality. In Figure 7.1, I plotted the bifurcation function a* for a = 2 and A = 1.5. I decide to let 1/A2 = 0.01 so that A 2 = 100. Using Equation 7.2 for A 2 = 100, I find m = 10 118 cells/lOO^m 3. In calculating the value of m, I used the following conversion: x 1 0 1 5 [ ^ ] x 1 0 - 1 5 [ ^ ] x (6.02 x 1 0 2 3 ) [ ~ ] Lira6 L pg kg 1 molecule 1 1 r , 9 nM X Mwt [kD~o) X 6.02xl023molecules [ M x Ll x W ^ 109 nM . Mwt where Mwt is the molecular weight of either IL-1/3 or T N F - a . In the previous conversion I used: • 1 L = 10 1 5 Lim3. • 1 kg = 10 1 5 pg. • 1 kg = 6.02 x 10 2 3 kDa. • 1 mole = 6.02 x 10 2 3 molecules. • 1 mole = 1 M x 1 L. • 1 M = 109 n M . Substituting the value of m into Equation 7.2, I find A\ = 37.50. For consistency, I also derived the homogeneous steady state chemical concentrations ci = 206 nM and c 2 = 1100 nM. From Equation 7.5, I calculate L\ = 70 fim and L 2 = 140 \xm. In order to generate the numerical simulations shown in the previous chapter, I have to provide the domain length. I chose L = 700 \im from the biological literature (Itagaki et al. 1989). I also considered a second set of parameter values by increasing the number of microglia in the system so that m = 40 cells/lOO^m 3, keeping all the other dimensional parameter values unchanged. This gives Ax = 150 and A2 = 400. I summarize the dimensionless parameter values for System 6.8 in Table 7.11. 119 Parameter Set 1 Set 2 Ax 37.50 150 A2 100 400 ei 0.0367 0.0367 0.0367 0.0367 a 2 2 A 1.5 1.5 domain length 5 5 Table 7.11: Dimensionless parameter values for the chemoattraction-repulsion model Equations 6.8 derived from biological parameter values given in Table 7.10. 2.2 Comparison with Biology In this section, I compare the theoretical results derived from the analysis of the chemo-attraction-repulsion model given by Equations 6.8, 6.10 and 6.11, with data from the biological literature. Thus, from the linear stability analysis, I can predict that the spacing between the plaques that could form under certain conditions is roughly given by L/n, where L represents the dimensional domain length and n is the mode. I have set L = 700 irni. In Figure 7.2, I show computer generated stationary spatial distributions of microglia and chemicals, using the two different sets of parameter values estimated from the experimental literature (as given in Table 7.11). The number of peaks (3, and respectively 5) observed in Figure 7.2a and 7.2b, correlates with modes n = 6 and n = 10, respectively. Therefore, I predict that the distance between the plaques in the domain is approximately 115 \im for the parameter set 1, and 70 pm for the parameter set 2. In order to compare the numerical results with the real biology (as found in Itagaki et al. (1989)), I summarize all data in Table 7.12. I observe that most numerical results 120 (a) Parameter set 1 in Table 7.11. (b) Parameter set 2 in Table 7.11. Figure 7.2: Stationary spatial distributions of microglia and chemicals for the nondimensional chemoattraction-repulsion model, using random initial conditions. (I used a log scale for the y-axis.) are in fairly good agreement with the real biology. For example, the distance between the plaques that was calculated from the numerical simulations falls within the range of values found in the biological literature, and some of the values also correlate well with my theoretical prediction of 115 \xm for the parameter set 1 and 70 [im for the parameter set 2. Note that I defined the size of a plaque as the distance over which the microglia density and the chemical concentration decay by a factor of e _ 1 . This definition is, to some extent, arbitrary, since there is no sharp boundary between what represents a "plaque" and what represents "the distance between the plaques" in Figure 7.2. Aspects to be Numerical Results Biology compared Set 1 Set 2 Plaques per 700 \xm 3 5 0-10 Sizes of plaques 80, 105, 140 pm 40, 70, 84, 105 fim 10-100 Lim Distance between plaques 105, 140 fim 40, 55, 105, 115 \im 50-200 Lim Table 7.12: Comparison between theoretical results and data from the biological literature (from Itagaki et al. (1989)). 121 54.6 [-0 I 1 1 1 1 1 0 1 2 3 4 5 space (b) Parameter set 2 in Table 7.11. Figure 7.2: Stationary spatial distributions of microglia and chemicals for the nondimensional chemoattraction-repulsion model, using random initial conditions. (I used a log scale for the y-axis.) 121a Additional biological interpretations can be derived by analyzing Figure 7.2. The differ-ence between the two parameter sets shown in Table 7.12 stems from the different values of the dimensional microglia homogeneous steady state, i.e. I took m = 10 cells/100/im 3 in set 1 (Figure 7.2a) and m = 40 cells/100/im 3 in set 2 (Figure 7.2b). From Itagaki et al. (1989) it is known that the number of reactive microglia cells correlates with the density of plaques, i.e. the higher the number of microglia cells, the more plaques develop. Thus, the numerical simulations are in good agreement with the biology. Also, microglia cells are concentrated at the center of the plaque. The chemicals have a higher concentration in the middle of the aggregate and they decay at the periphery. This type of cell and chemical distribution within a plaque is biologically supported by Itagaki et al. (1989). 2.3 Limitations of the Model Although the chemoattraction-repulsion model was successfully applied to mimic plaque development in Alzheimer's disease, I have to note many of its limitations. I summarize them as follows: • The model neglects several important biological aspects of Alzheimer's disease, such as the effects of neurons and other cell types (e.g. astrocytes). Many different chemicals secreted or uptaken by these cells, as well as various chemicals not even described in Chapter 2, were also omitted. • The model describes only one functional state of microglia cells, in particular re-active microglia. The distinction between the three states, resting, reactive and phagocytic, should also be modeled. In addition, the model should account for cell division and death. • The microglia random motility coefficient / i , the chemoattraction and chemore-pulsion coefficients Xi a n d X2, and the chemical diffusion coefficient D, were all assumed constant. This assumption was due to neglected changes occurring in 122 tissue properties, changes that would lead to functional forms of these coefficients. The model is one dimensional rather than two or three dimensional, and therefore cannot reproduce the round form of plaques, as seen in Figure 2.3 for example. Numerical simulations generated for a two dimensional model could lead to the formation of spots, as seen in Myerscough et al. (1998). The model only describes certain stages of plaque formation, in particular the formation of diffuse plaques. Thus, I should also model the progression from diffuse plaques to dense plaques. Even though I have used random initial conditions to generate all the numerical simulations, the obtained distribution of aggregates is sometimes periodic (Fig-ure 7.2a). This does not correlate well with the partially-ordered distribution seen in Figure 2.3, for example. The periodicity is the result of an artifact of the features built into the chemoattraction-repulsion model. The model is quite restrictive in fitting parameters that lead to plaque formation. For example, a small change in the value of the decay rate b\ or b2 (which easily occurs in lab experiments) could push us in region 1 of the parameter space A a (Figure 6.5) where there is no pattern formation. 123 Chapter 8 Conclusion 1 Concluding Remarks In this thesis I analyzed several chemotactic models developed for pattern formation in biological systems. My research was motivated by finding that microglia chemotaxis is important in plaque formations in the brains of Alzheimer's disease patients. I started with the study of the Keller-Segel model, a simple setting which incorporated the linear kinetics function g(m, c) = am—be. Analysis of stationary solutions established the existence of non-homogeneous steady state distributions of microglia and toxic chemical in the system. For different parameter sets, the dispersion relation indicated when the system was stable and when it was unstable. In particular, when instability was initiated by periodic initial conditions, the number of peaks that could form in the system was predicted. Numerical simulations using random initial conditions showed that the number of peaks correlated with the value of the most unstable mode given by the system. I should note that generating numerical simulations was a challenge, and a lot of effort was put into this task. The encountered problem was due to the well-known formation of very sharp peaks of cell densities (or chemical concentration) in any chemotactic system. Initially, I used the X P P simulation package which was not able to follow these sharp gradients. I then successfully used an adaptive mesh program developed by R. Russell, Simon Fraser University. 124 In addition, the linear stability analysis of the Keller-Segel model led to an aggregation condition which suggested possible therapy interventions for Alzheimer's disease. More specifically, I found that decreasing the number of reactive microglia or the production rate of /3-amyloid could slow the progression of plaque formation in Alzheimer's disease. The analysis continued with the examination of more complex kinetics function. I looked at the dynamics of the Keller-Segel model when the limiting kinetics function g(m, c) = ac(l — c) — bc(cCT\t — c)e~ac was used in the chemical equation. When I constructed this function, I took into consideration that microglia cells should be able to either produce or remove the toxic chemical from the system, depending on the parameter values; I also wanted to limit the formation of sharp peaks of microglia and chemical. The results of the numerical simulations generated for this particular chemotactic model were unexpected. Thus, the time evolution of microglia spatial distribution showed that, at the beginning of the simulation, the number of clusters correlated with the most unstable mode n; later, only one or two peaks could be seen in the domain. I concluded that this behaviour was probably induced by the nonlinear terms present in the equations. In addition, I am still investigating the reasons behind the different stationary cell distributions developed in this model and the Keller-Segel model. I also studied the chemoattraction-repulsion model for chemotaxis of cells in response to a chemoattractant and a chemorepellent that included a greater level of detail. The model predicted that local attraction and long range repulsion create periodic pattern forma-tion. Biological interpretations of the linear instability condition led to suggestions how therapy interventions in Alzheimer's disease might affect plaque distribution or tendency of plaques to form. For example, an increase in the production rate of chemorepellent (TNF-a), or a decrease in the level of activated microglia, could improve the symptoms of the disease. Real, biological parameters were successfully used in the simulations of this model. Stationary distributions of microglia and chemicals indicated that the number of 125 plaques developing in the domain was given either by the most unstable mode (when the values of unstable modes were close to the bifurcation mode) or by an arbitrary mode (when the values of unstable modes were significantly higher than the bifurcation mode). A final chapter was dedicated to the derivation of the biological parameters used in the numerical simulations of the chemoattraction-repulsion model, and to make the link be-tween this particular model and the biology of Alzheimer's disease. The parameters were either directly found in the literature or were calculated utilizing data from the literature. When none of the above was possible (as in the case of the chemorepulsion coefficient X2): I estimated the missing parameter using the regimes predicted by the model for appropriate phenomena. Based on the properties of the chemicals found in the brains of Alzheimer's disease patients and described in Chapter 2, I decided that microglia, IL-1/3 and TNF-cv should respectively play the role of the cell, the chemoattractant and the chemorepellent in the chemoattraction-repulsion model. I also verified that assumptions and conditions imposed on the model were reasonably well satisfied by the biological parameter values. Due to the ranges of the parameter values, any of the following three situations could occur: no pattern formation, one or two aggregates forming or periodic peaks forming in the domain. For my modelling purposes, I chose a set of parameters that mimic plaque formation in Alzheimer's disease. 2 Future Work The work presented in this thesis can be extended by continuing both theoretical and experimental analysis of the chemotactic systems introduced. I suggest the following research pathways for future investigation: • A nonlinear analysis should be developed to account for a more accurate interpre-tation of the system's dynamics far from the homogeneous steady state. 126 A chemotaxis model of two glial cells (microglia and astrocytes) and several at-tractant and repellent chemicals is already under investigation. Thus, I would like to incorporate additional biological variables into the systems in order to describe more realistic models. Higher dimensions of the chemotactic system should be studied. The brain is a stratified three dimensional system. Minimally, a layer (2D) with some thickness should be considered. Extensive studies should be done to include relevant biological parameters into my models, especially as the models increase in complexity. The correlation between theoretical model results and experimental observations will become, hopefully, stronger. In Figure 8.1, I show some of the complex interactions between the cells that play a role in Alzheimer's disease. Using this information, additional biological aspects, such as effects of neurons and various chemicals, could be analyzed and incorporated in the models. I should also take into consideration the difference between diffuse and dense plaques, between resting, reactive and phagocytic microglia. I could develop an expanded model to include some of the biological details de-scribed in the previous paragraph. For example, L. Edelstein-Keshet and A. Spiros have proposed a model that accounts for the health of neurons, activated and phagocytic microglia, astrocytes, and various chemicals (personal communication). Lab experiments could be initiated to test several predictions: identify effects of a variety of chemical gradients on microglia (attraction or repulsion); determine the chemoattraction and chemorepulsion coefficients Xi a n d X2 experimentally, and verify if they satisfy conditions stated in the analysis; investigate whether any real chemorepulsion can be detected in Alzheimer's disease. 127 Fibrillar igure 8.1: A schematic diagram of how neurons, microglia and astrocytes square symbols) interact with proteins (in round symbols) and tokines (in triangular symbols). Cel ls Cytokines M r Resting microglia IL-1/3 Interleukin-1/3 M a Active microglia IL-6 Interleukin-6 MPH Phagocytic microglia T N F - a Tumor necrosis factor-a Ay Resting astrocyte A s Stressed astrocyte Proteins N r Resting neuron ApoE Apolipoprotein E N s Stressed neuron A P P Amyloid precursor protein N d Dead neuron PAP yS-amyloid protein 128 Bibliography Ard, M . , G. Cole, J. Wei, A . Mehrle, and J. Fratkin (1996). Scavenging of alzheimer's amyloid /3-protein by microglia in culture. J o u r n a l of N e u r o s c i e n c e Research 43, 190-202. Bagnard, D., N . Thomasset, M . Lohrum, A. W. Piischel, and J. Bolz (2000, February). Spatial distributions of guidance molecules regulate chemorepulsion and chemoat-traction of growth cones. J o u r n a l of N e u r o s c i e n c e 20(3), 1030-1035. Banati, R. B. and K . Beyreuther (1995). Alzheimer's disease. In H. Kettenmann and B. R. Ransom (Eds.), Neuroglia;-Chapter 68, pp. 1027-1043. New York: Oxford University Press. Banati, R. B. , J. Gehrmann, P. Schubert, and G. W. Kreutzberg (1993). Cytotoxicity of microglia. G l i a 7, 111-118. Barron, K . D. (1995). The microglial cell, a historical review. J o u r n a l of the Neurolog-ical Sciences 134(Supp\.), 57-68. Benjamin, D., S. Wormsley, and S. Dower (1990). Heterogeneity in interleukin (il)-l receptors expressed on human b cell lines, differences in the molecular properties of il-1 alpha and il-1 beta binding sites. J o u r n a l of B i o l o g i c a l Chemistry 265(17), 9943-9951. Benveniste, E . N . (1995). Cytokine production. In H. Kettenmann and B. R. Ransom (Eds.), N e u r o g l i a , Chapter 5, pp. 700-713. New York: Oxford University Press. Boyce, W. E. and R. C. DiPrima (1992). Elementary Differential E q u a t i o n s and B o u n d -ary Value Problems (Fifth ed.). New York: John Wiley & Sons. Burden, R. L. and J. Faires (1989). N u m e r i c a l A n a l y s i s (Fourth ed.). Boston: PWS-K E N T Publishing Company. Chicoine, M . R., C. L. Madsen, and D. L. Silbergeld (1995). Modification of human glioma locomotion in vitro by cytokines E G F , bFGF, PDGFbb, N G F and T N F - a . Neurosurgery 36(6), 1165-1171. Corless, R., G. Gonnet, D. Hare, D. Jeffrey, and D. Knuth (1996). On the lambert w function. Advanced C o m p u t a t i o n a l Mathematics 5, 326-359. Cowley, G. (2000, January 31). Alzheimer's: Unlocking the Mistery. Newsweek, 46-51. Cross, M . and P. Hohenberg (1993). Pattern formation outside of equilibrium. Reviews of M o d e r n Physics 65(3), 851-10. Davis, J. B. , H . F. McMurray, and D. Schubert (1992). The amyloid beta-protein of alzheimer's disease is chemotactic for mononuclear phagocytes. B i o c h e m i c a l and B i o p h y s i c a l Research C o m m u n i c a t i o n s 189(2), 1096-1100. 129 de Castro, F., L. Hu, H. Drabkin, C. Sotelo, and A. Chedotal (1999). Chemoattraction and chemorepulsion of olfactory buln axon by different secreted semaphorins. The J o u r n a l of Neuro science 19(11), 4428-4436. Dickson, D. W. (April 1997). The Pathogenesis of Senile Plaques. J o u r n a l of Neu-ropathology and Experimental Neurology 56(4), 321-339. Dickson, D. W., S. C. Lee, M. -L . Zhao, A . Fishman, and S.-H. C. Yen (1995). The role of microglia in amyloid deposition in Alzheimer's disease. In P. F. Zatta and M . Nicolini (Eds.), N o n - N e u r o n a l Cells in Alzheimer's Disease, Chapter 10, pp. 108-127. Singapore: World Scientific. Ding, A. , E. Sanchez, S. Srimal, and C. Nathan (1989). Macrophages rapidly internalize their tumor necrosis factor receptors in response to bacterial lipopolysaccharide. J o u r n a l of B i o l o g i c a l Chemistry 264(7), 3924-3929. Duffy, P. E. (1983). Astrocytes: Normal, Reactive, and Neoplastic. New York: Raven Press. Edelstein-Keshet, L. (1988). Mathematical Models in Biology. New York.: The Random House. Faber-Elman, A. , V . Lavie, I. Schvartz, S. Shaltiel, and M . Schwartz (1995). Vit-ronectin overrides a negative effect of tnf-a on astrocyte migration. F A S E B Jour-nal 5(15), 1605-1613. Faber-Elman, A . , A . Solomon, J. Abraham, M . Marikovsky, and M . Schwartz (1996). Involvement of wound-associated factors in rat brain astrocyte migratory response to axonal injury: In vitro simulation. J o u r n a l of C l i n i c a l Investigation 97(1), 162-171. Farrell, B. , R. Daniele, and D. Lauffenburger (1990). Quantitative relationships be-tween single-cell and cell-population model parameters for chemosensory migration responses of alveolar macrophages to c5a. C e l l Motility and the Cytoskeleton 16(4), 279-293. Fiala, M . , L. Zhang, X . Gan, B. Sherry, D. Taub, M . C. Graves, S. Hama, D. Way, M . Weinand, M . Witte, D. Lorton, Y . - M . Kuo, and A. E. Roher (1998, July). Amyloid-/? induces chemokine secretion and monocyte migration across a human blood-brain barrier model. M o l e c u l a r M e d i c i n e 4(7), 480-489. Fil l i t , H. , W. Ding, L. Buee, J. Kalman, L. Altstiel, B. Lawlor, and G. Wolf-Klein (1991). Elevated circulating tumor necrosis factor levels in alzheimer's disease. Neu-roscience Letters 129, 318-320. FitzGerald, M . (1992). Neuroanatomy: Basic and c l i n i c a l (Second ed.). London: Bailliere Tindall. Gehrmann, J., Y . Matsumoto, and G. W. Kreutzberg (1995). Microglia: intrinsic immuneffector cell of the brain. B r a i n Rerseach Reviews 20, 269-287. 130 Gerald, C. F. and P. O. Wheatley (1989). Applied N u m e r i c a l Analysis (Fourth ed.). USA: Addison-Wesley Publishing Company. Gierer, A . and H. Meinhardt (1972). A theory of biological pattern formation. Kiber-netik 12, 30-39. Gilad, G. and V . Gilad (1995). Chemotaxis and accumulation of nerve growth factor by microglia and macrophages. J o u r n a l of Neuroscience Research 41, 594-602. Giulian, D. (1995). Microglia and neuronal dysfunction. In H. Kettenmann and B. R. Ransom (Eds.), N e u r o g l i a , Chapter 44, pp. 671-684. New York: Oxford University Press. Goedert, M . (1993). Tau protein and the neurofibrillary pathology of Alzzheimer's disease. Trends in Neuroscience 16(11), 460-465. Goodhill, G. J. (1997). Diffusion in axon guidance. E u r o p e a n J o u r n a l of Neuro-science 9, 1414-1421. Gradon, L. and A. Podgorski (1995). Displacement of alveolar macrophages in air space of human lung. M e d i c a l & B i o l o g i c a l E n g i n e e r i n g & Computing 33, 575-581. Griffin, W. S., J. G. Sheng, and R. E. Mrak (1997). Inflammatory pathways. In W. Wasco and R. E. Tanzi (Eds.), M o l e c u l a r Mechanisms of Dementia, Chapter 10, pp. 169-176. Humana Press: Humana Press. Griffin, W. S. T., J. G. Sheng, M . C. Royston, S. M . Gentleman, J. E. McKenzie, D. I. Graham, G. W. Roberts, and R. E. Mrak (1998). Glial-neuronal interactions in alzheimer's disease: The potential role of a cytokine cycle in disease progression. B r a i n Pathology 8, 65-72. Grindrod, P., J. Murray, and S. Sinha (1989). Steady-state spatial patterns in a cell-chemotaxis model. IMA J o u r n a l of Mathematics Applied in M e d i c i n e & Biology 6, 69-79. Hammacher, A. , R. Simpson, and E. Nice (1996). The interleukin-6 (il-6) partial an-tagonist (ql59e,tl62p)il-6 interacts with the il-6 receptor and gpl30 but fails to in-duce a stable hexameric receptor complex. J o u r n a l of B i o l o g i c a l Chemistry 271 (10), 5464-5473. Hardy, J. (1997). Amyloid, the presenilins and Alzheimer's disease. Trends in Neuro-science-20(4), 154-159. Harrison, L. G. (1993). Kinetic Theory of L i v i n g Pattern. New York: Cambridge Uni-versity Press. Itagaki, S., P. L. McGeer, H. Akiyama, S. Zhu, and D. Selkoe (1989). Relationship of microglia and astrocytes to amyloid deposits of alzheimer's disease. J o u r n a l of Neuroimmunology 24(3), 173-182. Jones, N . (June 2000). Soothing the inflammed brain. Scientific A m e r i c a n , 11-12. 131 Jordan, J., M . F. Galindo, R. J. Miller, C. A . Reardon, G. S. Getz, and M . J. Ladu (1998, January). Isoform-specific effect of apolipoprotein E on cell survival and beta-amyloid-induced toxicity in rat hippocampal pyramidal neuronal cultures. J o u r n a l of N euro science 18(1), 195-204. Kandel, E. R., J. H. Schwartz, and T. M . Jessell (1991). P r i n c i p l e s of N e u r a l Science (Third ed.). Connecticut: Appleton & Lange. Keller, E. F. and L. A . Segel (1970). Initiation of slime-mold aggregation viewed as an instability. J o u r n a l of T h e o r e t i c a l Biology 26, 399-415. Keller, E. F. and L. A . Segel (1971a). Model for chemotaxis. J o u r n a l of T h e o r e t i c a l Biology 30, 225-234. Keller, E. F. and L. A . Segel (1971b). Travelling bands of chemotactic bacteria: A theoretical analysis. J o u r n a l of T h e o r e t i c a l Biology 30, 235-248. Kornyei, Z., A . Czirok, T. Vicsek, and E. Madarasz (2000). Proliferative and migratory responses of astrocytes to in vitro injury. J o u r n a l of N euro science Research 61 (4), 421-429. Krul , E. and J. Tang (1992). Secretion of apolipoprotein e by an astrocytoma cell line. J o u r n a l of N e u r o s c i e n c e Research 32(2), 227-238. Lauffenburger, D. A . and C. Kennedy (1983). Localized bacterial infection in a dis-tributed model for tissue inflammation. J o u r n a l of M a t h e m a t i c a l Biology 16, 141— 163. Lee, C , M . Hoopes, J. Diehl, W. Gilliland, G. Huxel, E. v. Leaver, K . McCann, J. Umbanhowar, and A. Mogilner (2001). Non-local concepts and models in biology. J o u r n a l of T h e o r e t i c a l Biology 210(2), 201-219. Lee, S., W. Liu, D. Dickson, C. Brosnan, and J. Berman (1993). Cytokine production by human fetal microglia and astrocytes. Differential induction by lipopolysaccha-ride and IL-1. J o u r n a l of Immunology 150(7), 2659-2667. Lockhart, B. , C. Benicourt, J.-L. Junien, and A. Privat (November 1, 1994). Inhibitors of free radical formation fail to attenuate direct beta amyloid 25_35 peptide-mediated neurotoxicity in rat hippocampal cultures. J o u r n a l of N euros cience Research 39(A), 494-505. Mackenzie, I. R., C. Hao, and D. G. Munoz (1995). Role of microglia in senile plaque formation. N e u r o b i o l o g y of aging 16(5), 797-804. Maini, P., M . Myerscough, K.H.Winters, and J. Murray (1991). Bifurcating spatially heterogeneous solutions in a chemotaxis model for biological pattern generation. B u l l e t i n of M a t h e m a t i c a l Biology 53(h), 701-719. Mark, M . D., M . Lohrum, and A. W. Piischel (1997). Patterning neuronal connections by chemorepulsion: the semaphorins. Cell Tissue Research 290, 299-306. 132 Marzolo, M . P., R. von Bernhardi, G. Bu, and N. C. Inestrosa (2000, May 1). Expres-sion of a2-niacroglobulin receptor/low density lipoprotein receptor-related protein (LRP) in rat microglial cells. J o u r n a l of N e u r o s c i e n c e Research 60(3), 401-411. Mazel, T., Z. Simonova, and E. Sykova (1998). Diffusion heterogeneity and anisotripy in rat hippocampus. N e u r o R e p o r t 9(7), 1299-1304. McGeer, E . G. and P. L. McGeer (1999). Brain inflamation in alzheimer's disease and the therapeutic implications. C u r r e n t P h a r m a c e u t i c a l Design 5, 821-836. McGeer, P. L . and E. G. McGeer (1995). The inflammatory response system of brain: implications for theraphy of alzheimer and other neurodegenerative diseases. B r a i n Research Reviews 21, 195-218. McLean, C., R. Cherny, F. Fraser, S. Fuller, M . Smith, K . Beyreuther, A . Bush, and C. Masters (1999). Soluble pool of /3-amyloid as a determinant of severity of neurodegeneration in Alzheimer's disease. A n n a l s of N e u r o b i o l o g y ^(5(6), 860-866. Meinhardt, H . (1978). Models for the ontogenetic development of higher organisms. Reviews of Physiology, Biochemistry and P h a r m a c o l o g y 80, 47-104. Merrill, J. E . and E. N. Benveniste (1996). Cytokines-in inflammatory brain lesions: Helpful and harmful. Trends in N e u r o s c i e n c e 19(8), 331-338. Mesulam, M . - M . (November 1999). Neuroplasticity Failure in Alzheimer's Disease: Bridging the Gap between Plaques and Tangles. N e u r o n 24 , 521-529. Michishita, M . , Y . Yoshida, H. Uchino, and K . Nagata (1990). Induction of tumor necrosis factor-a and its receptors during differentiation in myeloid leukemic cells along the monocytic pathway, a possible regulatory mechanism for tnf-a produc-tion. J o u r n a l of B i o l o g i c a l Chemistry 265(lb), 8751-8759. Moghe, P. V. , R. D. Nelson, and R. T. Tranquillo (1995). Cytokine-stimulated chemo-taxis of human neutrophils in a 3-d conjoined fibrin gel assay. J o u r n a l of Immuno-logical Methods 180, 193-211. Mogi, M . , M . Harada, T. Kondo, P. Riederer, H. Inagaki, M . Minami, and T. Nagatsu (1994). Interleukin- 1 8, interleukin-5, epidermal growth factor and transforming growth factor-o; are elevated in the brain from parkinsonian patients. N e u r o s c i e n c e Letters 180, 147-150. Morse, P. and H. Feshbach (1953). Methods of T h e o r e t i c a l Physics. New York: McGraw-Hill. Mrak, R. E. , J. G. Sheng, and W. S. Griffin (2000). Glial cytokines in neurodegenerative conditions. In Patterson, Kordon, and Christen (Eds.), N e u r o - I m m u n e Interactions in N e u r o l o g i c and P s y c h i a t r i c D i s o r d e r s , pp. 9-17. Berin: Springer Verlag. Mrak, R. E. , J. G. Sheng, and W. S. T. Griffin (1995). Glial cytokines in alzheimer's disease: Review and pathogenic implications. H u m a n Pathology 26, 816-823. Murray, J. D. (1993). M a t h e m a t i c a l Biology (Second ed.). Berlin, Heilderberg: Springer-Verlag. 133 Myerscough, M . , P. Maini, and K . Painter (1998). Pattern formation in a generalized chemotactic model. B u l l e t i n of M a t h e m a t i c a l Biology 60(1), 1-26. Nash, J. M . (July 17, 2000). The new science of Alzheimer's. T i m e , 32-39. Nicholson, C. and E. Sykova (1998). Extracellular space structure revealed by diffusion analysis. Trends in N e u r o s c i e n c e 21(5), 207-215. Nilsson, L . , J. Rogers, and H. Potter (1998, Apri l 16). The essential role of inflamma-tion and induced gene expression in the pathogenic pathway of alzheimer's disease. F r o n t i e r s in B i o s c i e n c e 3, 436-446. Nolte, C , T.- Moller, T. Walter, and H. Kettenmann (1996). Complement 5a controls motility of murine microglial cells in vitro via activation of an inhibitory G-protein and the rearrangement of the actin cytoskeleton. N e u r o s c i e n c e Science 75(4), 1109-1120. Othmer, H. G. and A. Stevens (1997). Aggregation, blowup, and collapse: the ABC' s of taxis in reinforced random walks. SIAM J o u r n a l of A p p l i e d Mathematics 57(4), 1044-1081. Painter, K . , P. Maini, and H. Othmer (1999, May). Stripe formation in juvenile po-macanthus explained by a generalized Turing mechanism with chemotaxis. Proceed-ings of the N a t i o n a l Academy of Sciences in USA 96, 5549-5554. Painter, K . , P. Maini, and H. Othmer (2000). Development and applications of a model for cellular response to multiple chemotactic cues. J o u r n a l of M a t h e m a t i c a l Biology 41, 285-314. Pennica, D., V . Lam, N . Mize, R. Weber, M . Lewis, B. Fendly, M . Lipari, and D. Goed-del (1992). Biochemical properties of the 75-kda tumor necrosis factor receptor, characterization of ligand binding, internalization, and receptor phosphorylation. J o u r n a l of B i o l o g i c a l Chemistry 267(29), 21172-21178. Press, W. H. , B. P. Flannery, S. A . Teukolsky, and W. T. Vetterling (1989). N u m e r i c a l Recipes in C. Cambridge: Cambridge University Press. Prigogine, I. and R. Levefer (1968). Symmetry breaking instabilities in dissipative systems, i i . J o u r n a l of C h e m i c a l Physics 48, 1665-1700. Rivero, M . A. , R. T. Tranquillo, H. M . Buettner, and D. A . Lauffenburger (1989). transport models for chemotactic cell populations based on individual cell behavior. C h e m i c a l E n g i n e e r i n g Science ^ (12 ) , 2881-2897. Sawada, M . , N . Kondo, A . Suzumura, and T. Marunouchi (1989). Production of tumor necrosis factor-alpha by microglia and astrocytes in culture. B r a i n Research 491, 394-397. Schaaf, R. (1984). Global behaviour of solution branches for some neumann problems depending on one or several parameters. J. fur die Reine und Angewandte M a t h e -matik 346, 1-31. 134 Schaaf, R. (1985). Stationary solutions of chemotaxis systems. T r a n s a c t i o n s of the A m e r i c a n M a t h e m a t i c a l Society 292(2), 531-556. Segel, L. A . (1984). M o d e l i n g dynamic phenomena in m o l e c u l a r and c e l l u l a r biology. Cambridge: Cambridge University Press. Selkoe, D. J. (1991, November). Amyloid protein and alzheimer's disease. Scientific A m e r i c a n , 68-78. Selkoe, D. J. (1994). Cell biology of the amyloid /3-protein and the mechanism of Alzheimer's disease. A n n u a l Review of Cell Biology 10, 373-403. Shaffer, L. M . , M . D. Dority, R. Gupta-Bansal, R. C.A.Frederickson, S. G. Younkin, and K . R. Brunden (1995). Amyloid /3 protein (a/3) removal by neuroglial cells in culture. N e u r o b i o l o g y of A g i n g 16(5), 737-745. Sheng, J. G. , R. E. Mrak, and W. S. T. Griffin (1998). Enlarged and phagocytic, but not primed, interleukin-1 a immunoreactive microglia increase with age in normal human brain. Acta N e u r o p a t h o l o g i c a 95(3), 229-234. Smits, H. , L.A.Boven, C. Pereira, J. Verhoef, and H.S.L.M.Nottet (2000). Role of macrophage activation in the pathogenesis of Alzheimer's disease and human im-munodeficiency virus type i-associated dementia. E u r o p e a n J o u r n a l of C l i n i c a l In-vestigation 30(6), 526-535. Stalder, M . , A . Phinney, A . Probst, B. Sommer, M . Staufenbiel, and M . Jucker (1999). Association of microglia with amyloid plaques in brains of APP23 transgenic mice. A m e r i c a n J o u r n a l of Pathology 154(G), 1673-1684. Streit, W. J. (1995a). Microglia cells. In H. Kettenmann and B. R. Ransom (Eds.), N e u r o g l i a , Chapter 5, pp. 85-96. New York: Oxford University Press. Streit, W. J. and C. A . Kincaid-Colton (1995b). The brain's immune system. Scientific A m e r i c a n November, 56-61. Streit, W. J., S. A . Walter, and N . A . Pennell (1999). Reactive microgliosis. Progress in N e u r o b i o l o g y 57, 563-581. Strogatz, S. H. (1994). N o n l i n e a r D y n a m i c s and C h a o s . Reading, M A : Addison-Wesley. Sykova, E. (1997). The extracellular space in the ens: Its regulation, volume, and geometry in normal and pathological neuronal function. The N euro scientist 3(1), 28-41. Sykova, E., T. Mazel, and Z. Simonova (1998). Diffusion constraints and neuron-glia interaction during aging. E x p e r i m e n t a l Gerontology 33(7/8), 837-851. Tranquillo, R., S. Zigmond, and D. Lauffenburger (1988). Measurement of the chemo-taxis coefficient for human neutrophils in the under-agarose migration assay. Cell Motility and the Cytoskeleton 11, 1-15. Turing, A . M . (1952). The chemical basis of morphogenesis. P h i l o s o p h i c a l T r a n s a c t i o n s of the Royal Society of L o n d o n B237, 37-72. 135 Venters, H. D., R. Dantzer, and K . W. Kelley (2000). A new concept in neurode-generation: Tnfcv is a silencer of survival signals. Trends in N e u r o s c i e n c e 23(4), 175-180. Wilkinson, P. C. and W. S. Haston (1988). Chemotaxis: An overview. In G. D. Sabato (Ed.), Methods in E n z y m o l o g y , Volume 162, Chapter 1, pp. 3-16. San Diego: Aca-demic Press, Inc. Wisniewski, H . M . and J. Wegiel (1995). Non-Neuronal Cells Involved in ^-Amyloid Deposition in Alzheimer's Disease. In P. F. Zatta and M . Nicolini (Eds.), N o n -N e u r o n a l Cells in Alzheimer's Disease, Chapter 15, pp. 194-203. Singapore: World Scientific. Yamaguchi, M . , M . Michishita, K . Hirayoshi, K . Yasukawa, M . Okuma, and K . Na-gata (1992). Down-regulation of interleukin 6 receptors of mouse myelomono-cytic leukemic cells by leukemia inhibitory factor. J o u r n a l of B i o l o g i c a l C h e m -istry 267(31), 22035-22042. Yankner, B. A . (May 1996). Mechanisms of Neuronal Degeneration in Alzheimer's Disease. N e u r o n 16, 921-932. Yao, J., L. Harvath, D. L. Gilbert, and C. A . Colton (1990). Chemotaxis by a CNS macrophage, the microglia. J o u r n a l of N e u r o s c i e n c e Research 27, 36-42. Zauderer, E. (1989). P a r t i a l Differential E q u a t i o n s of A p p l i e d Mathematics (Second ed.). New York: Wiley-Interscience Publication. Zhou, H. F., L. H.-C. Lee, and R. D. Lund (1990). Timing and patterns of astrocyte migration form xenogeneic transplants of the cortex and corpus callosum. The J o u r n a l of Comparative Neurology 292(2), 320-330. Zwillinger, D. (1989). Handbook of Differential E q u a t i o n s . San Diego: Academic Press, Inc. 136 Appendix A Mathematical Background 1 Derivation of the Transient Solution The initial boundary value problem consists of the diffusion equation dc _ d2c dt dx2 S = ^ S , (A.1) the boundary conditions dc dx and the initial condition = 0 , (A.2) x=0,L c{x, 0) = w(x) , 0 < x < L . (A.3a) I seek solutions of Equation A . l of the form: c(x, t) = c(x) + cT(x, t) , (A.4) where c(x) is the steady state solution of Equation A . l and cT(x, t) is called the transient function. Using the method of separation of variables as described in Boyce and DiPrima (1992), I first calculate the transient function as cT(x,t) = X(x)T(t) , (A.5) where X(x) is some function of x and T(t) is some function of t. I substitute A.5 into Equation A . l to obtain 1 V _ X" D~f ~ ~Y ' 137 Since x and t are independent of each other, each side must be a fixed constant, say k; hence, I can write l_T_ _x^_ ~b~T ~ ~Y ~ or, the equivalent system, ( V -kDT =0 (A.6) X" -kX =0 System A.6 represents a set of two linear ordinary differential equations since each equa-tion depends on a single variable. Thus, I have reduced the second order partial differ-ential equation A . l to a set of two ordinary differential equations A.6. Now I can solve each of these two ordinary differential equations by carefully considering three distinct cases for the constant k, i.e. k = 0, k > 0, and k < 0. If k = 0 then T'(t) = 0 , and T(t) = constant , which implies that the solution c(x, t) might only depend on x, case which I am not interested in. If k > 0 then V T and T(t) = e kDt The equation for X is X" - kX = 0 , 138 with its associated characteristic equation r 2 - k = 0 . I obtain the following solution for X(x): X(x) = d e ^ x + C2e~^x . (A.7) Calculating C\ and C2 using the boundary conditions A.2 leads to k = 0 which contradicts k > 0. Thus I am left with the case k < 0. It is general practice to rename k = —A 2 , where A > 0. Then 'T(t) = e-x2Dt and X(x) = deiXx + C2e~iXx . The complex exponentials give rise to X(x) = Ci cos(Ax) + c 2 sin(—\x) , ' (A.8) after renaming the constants. Using the boundary condition at x = 0 I obtain c2 = 0, and from the boundary condition at x = L, I find —CiAsin(AL) = 0 . In order to get a non-trivial solution for my problem, I have to take C i ^ 0 and let sin(AL) = 0. This relation is satisfied if TL7T A = A„ = — , n = l , 2 , 3 . . . , (A.9) which leads to the solution Xn(x) = c n cos(A n x) , (A.10) 139 Tn(t) = e - x » D t . ( A . l l ) The values of A n given by Equation A.9 are called the eigenvalues of the eigenvalue problem described by A . l and A.2. The corresponding non-trivial solutions X n ( x ) are called the eigenfunctions. I have found an infinite number of transient functions c^(x,t) = cncos(\nx)e- x™ Dtx , n = 1,2,3... , (A.12) which represent the fundamental solutions of Equation A . l . From the principle of su-perposition, it is known that any linear combination of d^(x, t) also satisfies the partial differential equation and the boundary conditions. Consequently, I assume that oo cT(x,t) = ^ 2 c n c o s ( \ n x ) e ~ x " D t , n = 1,2,3... . (A.13) 71 = 1 The last step is to find the unknown coefficients cn by adding the fundamental Solutions A.12 such that the initial condition c(x, 0) = w(x) is satisfied. Substituting the sum A.13 into the initial condition gives w(x) = '^^cncos(—x) . (A-14) 71=1 If w(x) and its derivative are continuous (or piecewise so) on 0 < x < L, then on this interval w(x) can be represented by an infinite series of cosines (Boyce and DiPrima 1992). I multiply Equation A.14 by cos( n^-x), integrate from 0 to L and use the orthogonal property of the cosine function on the right hand side of Equation A. 14, i.e. / cos (— x ) cos( —— x) ax Jo L L After integrating from 0 to L, I obtain 0 , m ^ n ^ , m = n L j w(x) cos(^-x) dx = cn 2 , 140 which gives the form of the coefficients cn as ' 0 This last relation concludes my search for the transient function cr{x,t). 2 fL nix zn = — \ w(x) cos(—x) dx , n = 1,2,3.. . L Jo L 2 The Steady State Solution I have to calculate the steady state c(x) of Equation A . l . For this, I let dt ' and integrate dx2 twice to obtain c(x) = CiX + c 2 . The boundary condition A.2 at x = 0 gives C\ = 0. c 2 cannot be exactly determined given my initial conditions and boundary conditions. For now, let c 2 = c and thus the steady state is homogeneous and it is given by c(x) - c . (A.15) 3 Derivation of the Finite Difference Equations I have the initial boundary value problem given by the partial differential equations dm d2m d ( dc\ . . ^„ . ST = " a y ^ & r & J • ( A - 1 6 A ) | = ^ + * K 0 , (A.16b) 141 with boundary conditions dm dx dc dx x=0,L x-0,L = o, o, and initial conditions m(x, 0) = u(x) , c(x, 0) = w(x) , 0 < x < L , 0 < x < L . (A.17a) (A.17b) (A.18a) (A.18b) In general, I assume that the solution function f(x) has a continuous fourth derivative. I expand the function f(x) in a third order Taylor polynomial about the point xn and evaluate at xn + h and xn — h: f(xn + h) = f(xn) + f(xn)h + J-^h2 + f"(Xn)u2 , f " ' M h 3 + f I V ^ 1 ) 6 24 xn < Ex < xn + h , -h4 f(xn - h) = f(xn) - f'(xn)h + f " ( x n ) u 2 f " ' ( x n ) u 3 , 6 24 Xn ~ h < & < Xn . -h4 (A.19) (A.20) Subtracting these two equations I obtain f(xn + h) - f(xn - h) =f(xn) + 2h J v " 7 6 xn — h < £, < xn + h . Using the common notation for numerical approximations, I obtain fn+X fn—X f = J n 2h + 0(h2), (A.21) where f'n = f ' ( x n ) , f„ = f " ( x n ) , and /„ = f(xn) is given by the initial conditions. Similarly, I can derive the following formula for the second derivative: f(xn + h)- 2f(xn) + f(xn - h) = f"(xn) + 12 h2 xn — h < 3, < xn + h , 142 or, using numerical approximations, f „ = fn+l - 2 / n + / „ _ ! + Q ^ ( A Equations A.21 and A.22 are called central difference formulas for the first and second order derivatives. When the boundary conditions involve derivatives, the domain must be extended beyond the interval, using fictitious points outside the interval (Gerald and Wheatley 1989). I let the interval [0, L] be divided into 50 equal subintervals to have n = 51 discrete points. Each variable in System A.16 is unknown at the points m(xi) = and C(XJ) = c*, for i = 0 , . . . , 50. I have to extend the domain [0, L] to include one extra point at each end to incorporate the no flux boundary conditions. Therefore, the unknown points are now m(xi) = rrii and c(xi) = Ci, for i = — 1 , . . . , 51. I can use Formulas A.21 and A.22 to replace the equations in System A.16 and the boundary conditions A. 17. The cell boundary condition A. 17a can be replaced by , mi - m_i m » = -^h— = °' , m 5 i - m 4 9 ™ 5 0 = - ^ ^ = 0 , to obtain mi = m_i , (A.23a) m5i = m 4 9 . (A.23b) ci = c_i , (A.24a) c 5 i = c 4 9 • (A.24b) Calculations using A.23 and A.24 lead to ^ ( 0 ) = 2 ^ ^ ^ - x ( m i + m ; ) 2 ( C l ~ C ° ) , (A.25a) 143 Similarly, I get dm m i + i - 2mj + m ^ i (mi+i + mj ) ( c i + i - Cj) - (m^ + m i _ 1 )(c i - Cj_i) i = l , . . . , 4 9 , (A.25b) ^ m f c : n \ o m 4 9 ~ m 5 0 ( ^ 4 9 + m 5 0 ) ( c 4 9 - C 5 0 ) , . \ ^ ( 5 0 ) = 2 „ - ^ X . (A.25c) Similarly, the chemical equation can be expressed as follows: ^ ( 0 ) = 2D^j^ + g(m0,c0) (A.26a) = D ^ - ^ + ^ + g(mi, a), i = 1,..., 49 , (A.26b) ^ ( 5 0 ) = 2DC-^^ + g(m50,c50) . (A.26c) The set of Equations A.25 and A.26 constitutes the full discretization of System A.16. 144 Appendix B The Keller-Segel Chemotaxis Mode l 1 The Jacobian Matrix Let us have a general two dimensional nonlinear system of ordinary differential equations |T = F(x,y), (B.la) ^ = G{x,y), (B.lb) where F(x, y) and G(x, y) are nonlinear functions. I assume that System B . l has a steady state (x,y), i.e. it satisfies F(x,y) = G(x,y) = 0. ' (B.2) Using the method of smal l perturbations (Edelstein-Keshet 1988), I linearize System B . l . Thus, let x' and y' be small perturbations of the steady state such that x{t) = x + x'(t) , (B.3a) y(t) = y + y'(t). (B.3b) Substitution of B.3 into System B . l leads to °-g = F(x + x',y + y'), (B.4a) ^ = G(x + x',y + y'). (B.4b) 145 I expand the nonlinear functions F and G in a Taylor series about the point (x, y). After dropping the prime signs, I obtain doc — = F(x,y) + Fx(x,y)x + Fy(x,y)y + 0(x2,y2,xy) , (B.5a) ^ = G(x,y) + Gx(x,y)x + Gv(x,y)y + 0(x2,y2,xy) . (B.5b) Keeping only the linear terms and neglecting the higher order ones, and using B.2, I find the linear system dx dt dy Fx(x,y)x + Fy(x,y)y , .(B.6a) = Gx(x,y)x + Gy(x,y)y . (B.6b) I associate with the linear System B.6 the following matrix of coefficients j _ ( an an \ _ ( Fx(x,y) Fy(x,y) \ \ a 2 i a 2 2 / \ Gx(x,y) Gy(x,y) ) J is called the Jacobian of the System B . l . Stability of the steady state (x, y) can be determined from the Jacobian J after I define the following quantities r = trace J = an + a 2 2 , (B.8) A = det J = a n a 2 2 - a12a2i . (B.9) The steady state (x, y) can be classified as follows (Edelstein-Keshet 1988; Strogatz 1994; Murray 1993): 1. saddle point for A < 0. 2. unstable node for A > 0, r > 0 and r 2 > 4A. 3. star nodes and degenerate nodes for r 2 = 4A (they can be stable or unstable). 146 4. unstable spiral for r > 0 and r 2 < 4A. 5. neutral stable center for A > 0 and r = 0. 6. stable spiral for r < 0 and r 2 < 4A. 7. stable node for A < 0, T < 0 and r 2 > 4A. 2 The Lambert W Function The transcendental equation EV = -4-c* (B.10) can be solved using Lambert W function (Corless et al. 1996). By definition, the solution of yey = v, ( B . l l ) where v is a constant, is y = W(v). (B.12) The function W(v) is called the Lambert W function. If v is a real number, then for —1/e < v < 0 there are two possible real values for W(v) (Corless et al. 1996). Both values can be easily found, for a given v, using Maple. One value satisfies W(v) > —1 and is found using LambertW(v); the other value satisfies W(y) < —1 and is calculated using LambertW(—1, v). Mathematical manipulation of Equation B.10 leads to X * -*c- XaM c e ^ = — . (B.13) / i pb If I let - - c * , (B.14) 147 (B.15) then c* is given by c = J/ , X where y = W(y) and is determined using Maple. 148 Appendix C Extensions of the Keller-Segel Model 1 Different Kinetics Functions for the Chemical Equa-tion The model equations are (C.la) + g(m,c) , (C.lb) where g(rn,c) = s(c) — mh(c). 1.1 Reduced Systems The space independent system is: dm ~dt dc dt g{m,c) . 0 (C.2b) (C.2a) The time independent system is: (C.3a) 0 = D + g{m,c) . (C.3b) 149 1.2 Linear Stability Analysis of the Space Free System The homogeneous steady state (m,c) of the System C l satisfies the space independent System C.2, and therefore, it follows that 9{m,c) = 0, (CA) or, its equivalent form, m = — - , (C.5) h(c) provided h(c) ^ 0. Stability of the homogeneous steady state (m, c) depends on the Jacobian J of the System C.2. The Jacobian J is / 0 0 \ J=l , (C.6) \ 9m 9c J where gm = j & ^ y 9c = fc\{jn>5y I define r = trace J = gc , (C-7) A = d e t J = 0 . (C.8) Since A = 0, the homogeneous steady state (m, c) is a stable node if r < 0, i.e. gc < 0. 1.3 Linear Stability Analysis of the Full System Stability of this spatially uniform steady state is determined using the method of small perturbations (Edelstein-Keshet 1988). In a neighbourhood of the homogeneous steady state (m, c), I let m(x,t) = fh + m'(x,t), c(x, t) = c + c'(x, t) , where rn'(x,t) and c'(x,t) are small perturbations from the steady state. 150 System C l becomes dm! d2m' fdm'dc' , , . d 2 c ' \ , _ n , at -a* H ^ + ( m + m ) ^ J ' ( a 9 a ) dc' d2c' ^ == D ^ + g(m + m',c + c'). (C.9b) I linearize System C.9, as shown in Appendix B l , to obtain, after I drop the prime signs, dm a2m _ d2c ,~ ^ „ . -al = ^ - ^ 6 V ' ( a i 0 a ) dc a2 c — = D — + gmm + gcc. (C.lOb) I look for solutions of the form m(x,t) = M e " cos(gx) , (C. l la) c{x,t) = Ceatcos(qx) . (C . l lb ) Substituting C . l l into System C.10, I get a system of algebraic equations in M and C: (a + i i q 2 ) M - x m q 2 C = 0 , (C.12a) - g m M + (a + D q 2 ~ g c ) C = 0 , (C.12b) which has a non-trivial solution if ( a + uq2 - y m g 2 \ = 0 . (C.13) -gm a + D q 2 - gc ) This leads to the dispersion relation a2 + B a + 7 = 0 , (C.14) where 0 = W 2 + Dq2-gc, (C.15) 7 = q2[ii{Dq2 - gc) - xrhgm] • (C.16) 151 The aggregation condition holds whenever a > 0 which leads to X m g m > p{Dq2 - gc) . (C.17) I calculate the partial derivatives of g(m,c) at the homogeneous steady state (m,c): 9m = ~h(c) , gc = s'(c) - mh'(c) . (C.18a) (C.18b) Inequality C.17 can be rewritten using C.5 and the derivatives C.18, as follows: S ( C ) / L . / - W D q 2 - s'(c) + ^$-h'(c) h(c) which after simplification and division by x leads to -s (c > ^ X L D q 1 - s'(c) + s(c) s NJi'(c) /i(c) Rearrangement of this last relation gives \_h,{c) fi] The left-hand side (LHS) of this inequality is equivalent to: ~X , be~ac{ccrit -2c- accent + cue*2) LHS = a( l - 2c) - ac(l - c) L A * &c(c c r i t - c)e" a(l — 2c) — ac(l — c) a( l — c) ^ X —\ (Ccrit - c) + (-c - ojcc c r i t + ac*2) Ccrit C -1 - ac c r i t + ac) = a( l — 2c) — ac(l — c) a( l — c) — ac(l — c) A* ccrit — c X 1 — c = — ac — ac(l — c)—h ac + aac(l — c) A* cc rit c = ac 1^ c cr i t ) /.. s / X - IT- + (1 - c) a [Ccrit - C) \ A* (C.19) (C.20) 152 2 Multiple Kinetics Function for the Chemotactic Model dm d2m d ( dc\ . ._, n i . = M-^v-X^- m — + / ( m , c ) , (C.21a) The model equations are 9c _ Ddr21 dt dx 2.1 Reduced Systems The space independent system is: d G D^+g(m,c). (C.21b) 9m ~dt dc di = f(m,c), (C.22a) = <?(m,c). (C.22b) The time independent system is: 0 = ^-4x{mdi)+f^ (a23a) 0 = D ^ + g(m,c). (C.23b) 2.2 Linear S tab i l i ty Ana lys i s of the Space Independent Sys tem The homogeneous steady state (m, c) of the System C.21 satisfies the space independent System C.22, and therefore, it follows that f(m,c) = 0 , (C.24a) g(m,c) = 0 . (C.24b) Stability of the homogeneous steady state (m, c) depends on the Jacobian J of the System C.2. The Jacobian J is J = | f m f c | , (C.25) 9m 9c 153 where fm = - , fc = if \,- -s, Qm = %?-\,- -M and qc = I s I , . . I define "' ^m\(m,c), J c oc\(rn,c), •3"1 am\(m,cy a L ac\[m,c) r = trace J = fm + gc, ' (C.26) A = det J = fmgc - fcgm . (C.27) The homogeneous steady state (rri, c) is stable provided r < 0 and A > 0 are simultane-ously satisfied, that is, if fm + gc < 0 , (C.28) fmgc-fcgm > o . (C.29) 2.3 Linear Stability Analysis of the Full System Stability of this spatially uniform steady state is determined using the method of small perturbations (Edelstein-Keshet 1988). In a neighbourhood of the homogeneous steady state (ffi, c), I let m(x,t) = ffi + m'(x, t) , c(x, t) = c + c'(x, t) , where m'(x,t) and c'(x,t) are small perturbations from the steady state. System C.21 becomes dm' d2m' d dt ^ dx2 Xdx (m'tO + f(m + m',c + c') , (C.30a) dc1 d2c' °^ = D ^ + g(m + m',c + c') . (C.30b) I linearize System C.30, as shown in Appendix B l , to obtain, after I drop the prime signs, din d^Tn d^c — — ~dt = V ^ - X f h - ^ + f m m + fcC, (C.31a) dc n^ 2° di dx'2 = D ~ + gmm + gcc. (C.31b) 154 I look for solutions of the form m(x,t) = M e a t cos(qx) , ' (C.32a) c(x,t) = Ceatcos(qx). (C.32b) Substituting C.32 into System C.31, I get a system of algebraic equations in M and C (a + M 2 - f m ) M - ( X m q 2 + f c ) C = 0 , (C.33a) - g m M + (a + D q 2 - g c ) C = 0 , (C.33b) which has a non-trivial solution if ( a + uq2 — fm —Yfhq2 — fc \ = 0 . (C.34) -gm a + D q 2 - gc J This leads to the dispersion relation CT2 + /?CT + 7 = 0 , (C.35) where 8 = nq2 + Dq2-gc, (C.36) 7 = L i D q 4 - (xmgm + Dfm + pgc) q2 + f m g c - f c g m . (C.37) 155 Appendix D The Chemoattraction-Repulsion Mode l 1 Derivation of the Nondimensional System The system to be nondimensionalized is: dm d2m d ( dc\ dc2 \ _ N -m = ^ - Yx ( , X i m ^ " X 2 m ^ ) ' ( D - l a ) dc • d2 c = Di-^-j + atm - biCi i = l,2. (D.lb) I introduce the following nondimensional quantities x t m Ci x — — , i — ^ , m — , Cj — , x t m Ci where x, i, m, and c\ are scales to be defined later. Using the newly introduced nondi-mensional parameters, System D . l becomes 9m* = f r i & i r f __?_(Xi^m*^I _ X 2 C 2 ^ dc*2\ ^ dt* x2 dx*2 dx* \ x2 dx* x2 dx* J ' dc* Dt d2c* Let L 2 =V^ D2762", L x = y / D i / h , r = L \ l l x . (D.3) Then, I define the scales x, i, m, and Ci as follows: x = L 2 , i = r, ci = ^1m, i = l,2. (D.4) bi 156 Because of the relationship between rh and c\, it follows that m = m , ci = ci , c2 = c 2 , (D.5) where (m, 61,62) represents the homogeneous steady state of the original System D . l . I substitute the Scales D.4 into System D.2, and I drop the asterisks to obtain the following nondimensional system: £2 dm ~dt dci ~dt dc2 dt d2m d dx2 dx d2ci Ai dci dx 4u dx dx2 d2C2 dx2 + a2(m — C i ) , + m - c 2 , where Ai = fib A , Di ' £2 = Do 2 The Method of Green's Function I apply the method of Green's function to solve equations 0 = d2^ dx2 ctiCi + Oi{m , i = l,2, a = L i (D.6a) (D.6b) (D.6c) (D.7) (D.8) where OL\ = a2 and a 2 = 1. The method of Green's function introduces a linear operator S such that Sci(x) = fi(x) , (D.9) where I define the following: S fi(x) dx2 aim(x) , a. (D.10) (D . l l ) 157 so that Equations D.8 and D.9 are equivalent. If I can find a function Gi(x,x') that satisfies . SGi(x,x') =6{x-x') , (D.12) then the desired solution can be expressed as /oo G i f a x ^ f i i x 1 ) ,dx' . (D.13) •oo Here S(x — x') is the Dirac delta function and x' is a parameter. Equation D.12 states that Gi(x,x') satisfies SGi(x,x') = 0 for all x except for x = x'. At x = x', I have an infinite spike concentration. Thus, I have to find Gi(x,x'). I first solve Equation D.12 in the two regions x < x' and x > x' where SGi(x,x') = 0, (D.14) The linear equation D.14 has exponential solutions Gi(x,x') = e r x , where r is calculated from the characteristic equation - (r 2 - a,) = 0 , (D.15) i.e. r = ±y/a~i. Equation D.15 yields Gi(x, x') = U i e ^ x + Vie~^x , for - oo < x < x' , (D.16a) and Gi(x, x') = U2e^x + V2e~^x , for x' < x < oo , (D.16b) where U{ and Vi are constants. To determine all the coefficients Ui and Vi, I have to use several properties of the function Gi(x, x'). Using the boundary conditions dCi/dx\0yL = 0, I will assume that when x —> ±00, dci/dx —>• 0. Since Gi(x,x') represents chemical 158 concentration, it follows that in region 1, where —co < x < x', approaches zero as x —» —co. Therefore, in region 1 I must have Vi = 0. Similarly, in region 2, where x' < x < co, G[(x,x') —» 0 as x —> co, which implies U2 = 0. So far I have Gi{x,x') U i e V 5 i x ^ for — co < x < x' V2e~^x , for x' < x < co . (D.17) To find Ui and V 2 , I use two conditions that must be satisfied by the Green's function D.17 at x = x'. First, Gi(x,x') is continuous at the point x', so that U i e ^ * * ' = V 2 e ~ ^ x ' (D.18) The second condition to be satisfied by Gi(x,x') is the so called jump condition at x'. Integration of Equation D.12 over the interval (—co, co) gives /oo /*oo S G i ( x , x') dx = / S(x — x') dx = 1 , -co J —oo (D.19) by properties of the Dirac delta function. Further substitution of D.10 into the left hand side of D.19 leads to G'i(x,x')\x=xi+ - G'i(x,x')\x=x>- = - 1 , (D.20) provided Gi(x,x') is continuous at x = x'. Equation D.20 represents the "jump condi-tion". The derivative of Gi(x,x') with respect to x is Gi{x,x') Ui y/c7ie^x , for — co < x < x' — V 2 J o 7 i e ~ ^ ° r i X , for x' < x < co . (D.21) Using Equation D.21, Equation D.20 becomes - V 2 J c 7 i e - ^ x ' - U x J b J i e ^ ' = - 1 . (D.22) 159 Equations D.18 and D.22 form a system of two equations with two unknowns U\ and V2 that can be easily solved, yielding 1 V2 1 2^/al -y/alx' Thus, I find that Green's function is given by i = e ^aix' e^aix for — oo < x < a:' i==ey/oTix e-^cTix for x' < a; < oo , (D.23) or, its equivalent form, Gi(x,x') 2,/cTi -^fai\x-x'\ (D.24) Finally, using Equations D . l l , D.13 and D.24, the chemical solution Q of Equation D.8 reads Ci(x) 'a. /oo e - y / « \ * - * ' \ m r x ' ) , dx' -oo (D.25) 3 The Shadow System In this section I derive the shadow system. I first rewrite Equation D.25 as d(x) = eJcTi{x x ' ) m . x . rOO ')dx'+ / e-^i{x-x,)m(x')dx' Jx' From D.26, I can calculate the derivative of Ci(x) with respect to x: dci _ y/oTj ~dx " ~Y~ a. y/a~ / e^x-x,)m(x') dx' - yfa~ e-^l{x-x')m(x') dx' — j sign(x' - x)e-^x-x'lm{x') dx1 . 2 (D.26) (D.27) 160 Substitution of Equation D.27 into the cell Equation D.6a leads to dm d2m d dt dx2 dx d2m d ( A i — f sign(x' - x)e-^x-x'\m{x') dx' /OO \ sign(x' - x)e-^x-x'lm{x') dx'J -co / m dx2 dx A2 r , , sign(x — xt ^ A e - a \ x - x ' \ _ e - | x - x ' | j m(x,j dx> m If I let ^ A />0O v = -^(K*m) = ^ K{x- x')m(x') dx' , K{x) = sign(-x) ( A e ~ a W - e - M ) } A = ^ l ? l £ l a = hi I obtain the shadow system dm d2m d . . ~m = ~dx^ ' d x { v m ) • (D.28) (D.29) (D.30) (D.31) 4 Linear Stability Analysis of the Shadow System I substitute the perturbed steady state m(x,t) = 1 + m' into Equation D.31 to obtain dm' d2m' d dt dx2 dx Al{K*(l+;m')){l + m') Li (D.32) Then, I have: /CO K(x - x')dx' + K *rri = K * rri , -co where I used the property /CO K(x' - x) dx' = 0 , -co since K(x) is an odd function integrated over a symmetric interval. Equation D.32 becomes dm' dm' d dt dx2 dx A , (K*m')(l + rri) 161 Keeping only the linear terms, I get dm! d2m' d ( A 2 M 8 * * V 2 ( / f * m ' ' j - ( D ' 3 3 ) Now let m'(x,t) = M e a t cos(gx). I need the derivative of K * m' with respect to x. I have the following steps: /oo K(x - x')Meatcos(qx') dx' •oo /oo K ( x - x') cos(gx') dx' -oo /oo X(x ' ) cos[g(x - x')] dx' -oo /oo K(x')[cos(gx) cos(gx') + sin(gx) sin(gx')] dx' -oo /OO /"OO K(x ' ) cos(gx') rix' + M e a t sin(gx) / i f (x') sin(gx') dx' OO J —oo poo = 2 M e a t sin(gx) / K(x ' ) sin(gx') dx' Jo = - 2 M e a t sin(gx) J ^Ae~ax' - e~x'^ sin(gx') dx' = - 2 M e a t sin(gx) ^ ^ g 2 + a 2 g 2 + 1 (D.34) where in the third step I used the commutative law of convolution, i.e. for two functions / and g, f * g = g * f (Zauderer 1989). I also observe that the kernel K(x) is an odd function, sin(gx) is also an odd function, while cos(gx) is an even function. Therefore, on a symmetric interval such as (—oo, oo), the integral of the product K(x) cos(gx) is zero. Now, (K * m') = - 2 M q e ^ cos(gx) - y - f - ^ - . (D.35) ox \ q2 + a2 q2 + 1 / Substituting m' = M e a t cos(gx) and Equation D.35 into Equation D.33, yields M a e a t cos(gx) = - q 2 M e a t cos(gx) + q 2 A 2 M e a t cos(gx) ( _ A , - J . (D.36) \ q2 + a2 q2 + 1 / 162 This is equivalent to a = - g 2 + g 2 A 2 ( - I A ^ - ^ - T ) . (D.37) -,2 _|_ a 2 „Z _|_ i 163
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Application of chemotactic models to Alzheimer’s disease
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Application of chemotactic models to Alzheimer’s disease Luca, Magdalena 2002
pdf
Page Metadata
Item Metadata
Title | Application of chemotactic models to Alzheimer’s disease |
Creator |
Luca, Magdalena |
Date Issued | 2002 |
Description | This thesis presents the analysis of several chemotactic models developed for pattern formation in biological systems. The project was motivated by the complex interactions between glial cells, cytokines and proteins that cause neuronal death and eventual dementia in Alzheimer's disease. I investigate conditions that lead to periodic patterns and aggregation of glial cells and chemicals, as observed in the senile plaques of Alzheimer's disease. I have examined a hierarchy of one dimensional models, starting with a model in which there is a single chemoattractant inflammatory agent and the microglia cells. Thus, I constructed a minimal partial differential equations model similar to the Keller-Segel model for the chemotaxis of microglia and their involvement in secretion and uptake of various inflammatory factors. Stationary solutions of the system of equations are investigated using a method developed by Schaaf. The linear stability analysis of the system leads to the aggregation condition which can be biologically interpreted. As a result, I suggest possible methods that might affect plaque formation in Alzheimer's disease. Numerical simulations help corroborate theoretical results. A different form of the kinetics function which is part of the Keller-Segel model is suggested. Thus, the limiting kinetics function in the chemoattractant equation incorporates production or removal of the chemical factor by the cells, depending on the model parameters. Also, it limits the sharp increase in the cell density which is well known to occur in this type of a chemotactic model. Finally, a model for chemotaxis of cells in response to a chemoattractant and a chemorepellent is developed to include a greater level of detail. The system is reduced in complexity using quasi-steady state approximations, and explicit solutions for chemical diffusion are found using the Green's function method. The model predicts that local attraction and long range repulsion create periodic pattern formation. Biological interpretations of the linear instability condition lead to suggestions for therapy interventions in Alzheimer's disease. Using real biological parameters available in the literature, I compare and discuss the applicability of this particular model to the pattern formation observed in Alzheimer's disease. |
Extent | 11253112 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-09-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0080045 |
URI | http://hdl.handle.net/2429/13098 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2002-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2002-731995.pdf [ 10.73MB ]
- Metadata
- JSON: 831-1.0080045.json
- JSON-LD: 831-1.0080045-ld.json
- RDF/XML (Pretty): 831-1.0080045-rdf.xml
- RDF/JSON: 831-1.0080045-rdf.json
- Turtle: 831-1.0080045-turtle.txt
- N-Triples: 831-1.0080045-rdf-ntriples.txt
- Original Record: 831-1.0080045-source.json
- Full Text
- 831-1.0080045-fulltext.txt
- Citation
- 831-1.0080045.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080045/manifest