Popular in Course
Popular in Mathematics (M)
This 190 page Class Notes was uploaded by Miss Hillary Grady on Thursday October 29, 2015. The Class Notes belongs to MATH 490 at College of William and Mary taught by Junping Shi in Fall. Since its upload, it has received 16 views. For similar materials see /class/231141/math-490-college-of-william-and-mary in Mathematics (M) at College of William and Mary.
Reviews for Seminar
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 10/29/15
Chapter 5 Reaction diffusion systems and pattern formation 51 Reaction diffusion systems from biology ln ecological problems different species interact with each other and in chemical reactions different substances react and produce new substances System of differential equations are used to model these events Similar to scalar reaction diffusion equations reaction diffusion systems can be derived to model the spatial temporal phenomena Suppose that Ut z and Vt z are population density functions of two species or concentration of two chemicals then the reaction diffusion system can be written as 8U 82U E DUW fFU7V7 51 8V 82V 39 E w t W where DU and DV are the diffusion constants of U and V respectively and FU V and CU V are the growth and interacting functions For simplicity we only consider one dimensional spatial domains The domain can either be the whole space foo 00 when considering the spread or invasion of population in a new territory or a bounded interval say 0 L In the latter case appropriate boundary conditions need to be added to the system Typically we consider the Dirichlet boundary condition hostile exterior Ut 0 Ut L 0 Vt0 Vt L 0 52 or Neumann boundary condition no ux Umt 0 Umt L 0 Vmt 0 Vmt L 0 53 The following are a list of models which has been studied by biologists and chemists 83 84 CHAPTER 5 REACTION DIFFUSION SYSTEMS Two species spatial ecological models Predator prey Ut DUUM AU 7 BU V V DVVm 7 CV kBU V 54 Here U is a prey and V is a predator AU is the growth rate of U CV is the death rate of V k gt O and BU V is the interaction between U and V Typical forms of these functions are AU aU or aU 17 2 CV CV dU JV 55 B U V dUV 7 0r 1 5U Other ecological models are Competition Ut DUUm AU 7 BU V Vt DVVM CV 7 kBU V 56 and Mutualism Ut DUUm 7 AU BU V V DVVm 7 CV 7 kBU V 57 Chemical reaction models See more introduction in Murray s book Thomas mechanism Ut DUUm a 7 U 7 pRU V Vt DVVM 01 7 V 7 pRU V 58 where RUV Thomas model is based on a speci c reaction involving the e fU gUZ substrates oxygen and uric acid which reacts in the presence of the enzyme uricase U and V are the concentration of the uric acid and the oxygen respectively and they are supplied at constant rates a and 0 Both U and V degrde linearly proportional to their concentrations and both are used in the reaction at a rate of pRU V RU V is an empirical formula For xed V RU V is nearly linear in U for small U and it approached to zero as U 7 oo Schnakenberg reaction U DUUm a 7 bU CUZV Vt vam 7 d 7 CUZV 59 This equation is from an arti cial tri molecular chemical reaction proposed by Schnakenberg X 7 A B 7 Y 2X Y 71 3X 510 k 71 Gierer Meinhardt model 52 TURING INSTABILITY AND PATTERN FORMATION 85 U2 U DUUma7bUc7 Vt DVVmdU27eV 511 Chlorine dioxide iodine and malonic acid reaction CIMA reaction 4UV UV UtDUUmm Slta7U7Wgt7 VtDvabltU71U2gt 512 Gray Scott model U DUUm 7 UZV f17 U Vt DVVm UZV 7 f kV 513 All these models can be rewritten to a standard form 814 82u 7d l 825 3m2 91quot after a nondimensionalization procedure 52 Turing instability and pattern formation We consider the nondimensionalized system ut umfuv tgt0 m 6 01 U dvm guv tgt 0 m 6 01 7 515 71020 um 1 M7570 Mt 1 07 u0z ax U0z We suppose that 140120 is a constant equilibrium solution Le 100107110 07 and 90407110 0 516 In fact the constant solution 140120 is also an equilibrium solution of a system of ordinary differ ential equations g uv tgt 0 517 In 1952 Alan Turing published a paper The chemical basis of morphogenesis which is now regarded as the foundation of basic chemical theory or reaction diffusion theory of morphogenesis Turing suggested that under certain conditions chemicals can react and diffuse in such a way 86 CHAPTER 5 REACTION DIFFUSION SYSTEMS as to produce non constant equilibrium solutions which represent spatial patterns of chemical or morphogen concentration It is noticeable that Alan Turing was one of the greatest scientists in twentieth century He was responsible for much of the fundamental work that underlies the formal theory of computation general purpose computer and arti cial intelligence During World War II Turing was a key gure in the successful effort to break the Axis Enigma code used in all German U boats which is one of the most important contributions scientists made to the victory of the war Turing s idea is a simple but profound one He said that if in the absence of the diffusion considering ODE 517 u and v tend to a linearly stable uniform steady state then under certain conditions the uniform steady state can become unstable and spatial inhomogeneous patterns can evolve through bifurcations In another word an constant equilibrium can be asymptotically stable with respect to 517 but it is unstable with respect to 515 Therefore this constant equilibrium solution becomes unstable because of the diffusion which is called di usion driven instability Diffusion is usually considered a smoothing and stablising process that is why this was such a novel concept More information of Turing instability can be found in Murray s book Chapter 14 So now we look for the conditions for the Turing instability described above If uo U0 is stable with respect to 515 then the eigenvalues of Jacobian J f f 518 gu 91 at uo 120 must all have negative real part which is equivalent to TraceJ fu 91 lt 0 Determinant fugv 7 fvgu gt 0 519 Now suppose that uo U0 is unstable with respect to 517 Then u0v0 must also be unstable with respect to the linearized equation 1 m fuu07 U0 gt fuU07U0117 75gt 07 90 E 071 11 dwm MAW U0 gt AMUO vow tgt 07 90 E 071 520 1t70 t71 t70 17571 07 gt07 90 C90 10 90 90 or in matrix notation Ilt DkIIm JII 521 Where 514535 and D 1 2 522 From the separation of variables for single equation we expect that j and 1 are both separable and we try the following form 11t 5 Wig W gtcosk7rz Wktcosk7rm 523 52 TURING INSTABILITY AND PATTERN FORMATION 87 We shall try this form of solution and the general solution for 11 would be 00 t 5 Z Wkt cosk7m 524 k0 By substitution of this form into 521 we obtain that Wkt satis es a system of ordinary differ ential equations de T T If 520 is an unstable system then at least one of Wkt would go to in nity as t 7 00 26 one of the eigenvalues of matrix Mk J 7 52sz has positive real part The eigenvalues p1 and p2 of J 7 2sz satisfy J 7 k27r2DWk 525 22 7 Af5gu7 d1k27erM 5 26 fu9v71011972ik2W2dfugvdk47r4 07 39 TraceMk fu 1 972 7 d 1k27r27 5 27 DeterminantMk L291 7 fugu2 7 k27r2dfu 971 dk47r439 39 Since fu gv lt 0 then TraceMk lt 0 is always true So Mk cannot be a source or spiral source If Mk is unstable then it can only be a saddle Therefore we must have DeterminantMk lt 0 Remember that we only need one of Mk is unstable So let s take a look of the values of hk27r2 DeterminantMk for different k or 15271392 indeed Here we de ne the function hw 7 595 7 1295 7 mm gm dwz 528 When w 0 710 fugv 7 fugu2 gt 0 because we require fugv 7 fvgu gt 0 and when 152 is large 7 00 we also have lim hw 00 So we must hope that for some 0 lt 152 lt 00 217200 hk27r2 lt 0 Since h is a quadratic equation we must require min hw lt 0 It is easy to calculate 71 dfu l 91 2d that min hw is achieved at we This implicitly implies that w dfu g gt 0 529 otherwise hw gt 0 for all w 2 O The minimum value at we is d 2 hw0 A2 595 7 fugu 7 hwy 530 which must be negative so dfu l 912 7 4dltfu9v 7 fugu gt 0 531 Summarizing all conditions we got so far if 140120 is stable with respect to 515 but unstable with respect to 517 then 1 and J must satisfy in l 91 lt 07 fug39u 7 fvgu gt 07 2 532 dfu 91 gt 07 91 7 4610072911 7 fugu gt 88 CHAPTER 5 REACTION DIFFUSION SYSTEMS The conditions in 532 give the range of parameter at such that uo 120 is possibly unstable but only 532 can only assure that certain values of hw are negative but it cannot guarantee there is a value k27r2 such that hk27r2 lt 0 To achieve that we directly calculate the conditions for hk27r2 lt O From 528 we have W279 129 7 fuguA2 7 k27r2dfu gm dk47r4 lt 0 533 and if fu gt 0 and gv lt 0 then ADA 7 kz n39zgv d gt ham 7 k27r2 E dkA 534 where D fugv 7 fugu dkA is a family of curves with dkA de ned on kZn39Zfu 00 see example below For a pair of parameters A d if d gt dkA for some natural number k then uo 120 is unstable But ifd lt dkA for any h then it is stable Summarizing the above calculation we conclude Theorem 51 Suppose that uouo is a constant equilibrium solution of the reaction di usion system ut U Afuu t gt 0 m 6 01 1 allm Aguu t gt 0 m 6 01 535 we 7 um 1 7 mt 0 7 mt 1 7 0 M0790 495 110795 M907 and the Jacobian of fg at u0U0 is Ju0u0 lt 5 71 If J at and A satisfy U 39U fugt079ult07 fu9ult07 fugu fu9ugt07 536 din gt gt 0 din gt 7 4dltfugi 7 firm gt 0 39 and d gt dkA for some k where A DA 7 l 27r2 dk 91 D fugu 7 qum 53937 k27r2qu 7 k27r27 then u0U0 is an unstable equilibrium solution with respect to 535 but a stable equilibrium solution with respect to u Afu U 1 Aguu We use the example of Schnackenberg to illustrate the theorem Consider the following and a nondimensionalized system of reaction diffusion equations describing the chemical reaction is ut ummAa7uu2u tgt0 m 6 01 at alum Ab 7 7m 25gt 0 w 6 01 uwt0 umt 1 umt0 umt 1 0 u0x A0m u0z B0m 538 52 TURING INSTABILITY AND PATTERN FORMATION 89 where a b gt 0 There is a unique equilibrium constant solution 140120 such that b b 539 U0 1 7 U0 14 and the Jacobian at 140120 is b 7 a b a2 Ju0v0 7 5 251 540 7 71 12 b L Then from simple computation the conditions 536 become 0ltb7alt ba3 ba2gt0 lt gt lt gt 541 11 7 a gt b a3 db 7 17 b 13 gt 4db a4 For simplicity and transparency of the result we select a pair of a b 1 2 then the rst row of 541 is satis ed The formula for dk in this case is 27 7 kW dkw k27r2 7 3k27r239 542 0 200 400 x 600 800 Figure 51 Parameter space d for Turing instability The rst three dk are drawn in Fig 51 Each dk is a concave up curve with a unique minimum point From Theorem 51 if 1 falls below all of dk then 140120 is stable for that parameter choice and if d is above any of dk then 140120 is unstable When 140120 is unstable the solution of linearized equation 520 has a form 0 Ilt m 2clkemktW1k CgkeMktWZk COSk7rm 543 k 0 90 CHAPTER 5 REACTION DIFFUSION SYSTEMS where pm and ngk are the eigenvalues of AJikzsz and Wlk ng are corresponding eigenvectors Since it is unstable then one or more eigenvalues ulk or 2 are positive For example if d is a point such that d gt d2 but 1 lt dk for all other k then there is one eigenvalue 12 gt 0 or complex with positive real part and all other eigenvalues have negative real part and as t a 00 Ilt z z 0125 12tW12t cos27rm 544 So the solution grows exponentially as t a 00 and it grows in a temporal spatial pattern W1kt cos27rm which is called an unstable mode of this equilibrium solution For an unstable equilibrium there may be more than one unstable modes if the parameter pair d is below two curves dk The linear analysis of Turing instability is completed now For nonlinear systems if the initial value is near the constant equilibrium solution then in a short time the orbit would be close to that of linearized equation When the equilibrium is unstable a spatial pattern is formed following its unstable mode and the pattern grows exponentially fast for a short period But after that the nonlinear structure takes over and usually nonlinear functions force the solutions staying bounded like in the logistic model thus the solution approaches a bounded spatial pattern as t a 00 Although the mechanism of the nonlinear evolution is still not clear mathematically people believe that the spatial pattern formed in the earlier stage near equilibrium solution is kept which seems to be supported by numerical simulations Many problems are still topics of current research interests More biological discussion can be found in Murray s book page 388 414 53 Turing patterns in 2D and animal coat patterns 54 Chemotactic motion of microorganisms 55 Spatial epidemic model Chapter 5 Exercises 1 Consider the diffusive predator prey model 2 8 UDU8 UAU7BUV sgt0y 0L 59 3 DviCVDUX sgt0y 0L 545 U 8V 8V 39 ho EL 8740570 EL 07 WOW Uoy7 WOW VOW Use the substitution Math 490 Project presentations April 23 Friday Bifurcation diagrams of equilibrium solutions of nonlinear diffusion models Young He Lee Lena Sherbakov and Jacky Taber April 26 Monday Animal aggregation and nonlinear diffusion models Michael Deal Brian Van Hise and David Pluim April 28 Wednesday Numerical computations of reaction diffusion systems Zayd Khoury Chris Leonetti and Andrew Todd April 30 Friday Two patch diffusive fishery harvesting model Victoria Dyer Eileen Tschetter and David Weissenberger Modeling of Biological Systems A Workshop at the National Science Foundation March 14 and 15 1996 http www nsf govbiopubsmobsmobs htm Panel Members m Peter Kollman University of California San Francisco Co Chair Simon Levin Princeton University Alberto Apostolico University of Padova Marjorie Asmussen University of Georgia Bruce L Bush Merck Research Labs Carlos Castillo Chavez Cornell University Robert Eisenberg Rush Medical College Bard Ermentrout University of Pittsburgh Christopher Fields Santa Fe Institute John Guckenheimer Cornell University Alan Hastings University of California Davis Michael Hines Yale University Barry Honig Columbia University Lynn Jelinski Cornell University Nancy Kopell Boston University Don Ludwig University of British Columbia Terry Lybrand University of Washington George Oster University of California Berkeley Alan Perelson Los Alamos National labs Charles Peskin Courant Institute of Mathematical Sciences Greg Petsko Brandeis University John Rinzel National Institutes of Health Robert Silver Marine Biological Laboratory Sylvia Spengler Lawrence Berkeley Labs DeWitt Sumners Florida State University Carla Wofsy University of New Mexico The common theme of this report is the tremendous potential of mathematical and computational approaches in leading to fun damental insights and important practical benefits in research on biological systems Mathematical and computational approaches have long been appreciated in physics and in the last twenty years have played an ever increasing role in chemistry In our opinion they are just coming into their own in biology Our report also highlights computational issues that are com mon across biology from the molecular to the ecosystem Com puters are getting more powerful at a prodigious rate and in parallel the potential for computational methods to ever more complex systems is also increasing Thus it is essential that the next generation of biological scientists have a strong train ing in mathematics and computation from kindergarten through graduate school I MOLECULAR AND CELLULAR BIOLOGY A central organizing theme in Molecular and Cellular Biology is the relationship between structure of molecules and high level complexes of molecules and their function both in normal and aberrant biological contexts The connection between structure and function was most clearly illustrated in the paper that began Molecular Biology the elucidation of the structure of DNA by Watson and Crick 1 GENOME Since 1990 genomic sequence information has continued its exponential growth Sequencing technology is be ing applied directly to sequence diversity analysis and gene ex pression analysis via high throughput chip based automated as say systems This influx has changed both the questions that are asked as well as the range of the interactions considered 4 Significant work is required to develop data management systems to make these data not just retrievable but usable as input to computations and amenable to complex ad hoc queries across multiple data types Significant work is also required on tech niques for integrating data obtained for multiple observables at different scales with different uncertainties data fusion and for formulating meaningful queries against such heterogeneous data data mining For example it should be possible in the future to ask what differ ences to expect in the kinetic efficiencies of a signal transduction pathway across multiple individuals given the differences in the sequences of the proteins involved in the pathway Answering such queries will require improvements in data models hetero geneous database management systems multivariate correlation analysis molecular structure prediction constrained network mod eling and uncertainty management 2 PROTEIN STRUCTURE AND FUNCTION As the amount of genomic data grows three dimensional struc ture will provide an increasingly important means for exploiting and organizing this information Structure provides a unique yet largely unexplored vehicle for deducing gene function from sequence data Structure also links genomic information to bio logical assays and serves as a basis for rational development of bioactive compounds including drugs and vaccines Research opportunities in this area can be divided into four dis tinct categories experimental structure determination structure prediction structure exploitation of globular proteins and mod eling of membrane proteins where the determination of high resolution structures is much more difficult 3 SIMULATIONS Simulations of molecules of biological in terest use computational representations that range from sim ple lattice models to full quantum mechanical wave functions of nuclei and electrons If one has access to a macromolecular structure derived from NMR or X ray crystallography then one can begin with a full atom representation and fruitfully exam ine small changes in the system such as ligand binding or site specific mutation Again the goal is to reproduce and predict structure dynamics and thermodynamics Developments in both hardware and software for parallel com puting have played a major role However the longest time simulations that have been carried out are still 9 orders of mag nitude away from the typical time scale for experimental pro tein folding Simplified but realistic models for example using a 7 continuum treatment of the solvent Gilson et al 1995 could increase the time scale by 1 2 orders of magnitude Contin uum representations may more readily incorporated into Monte Carlo methods and thus allow large movements of the molecule during simulation Senderowitz et al 1996 In some cases the use of Langevin and Brownian dynamics and multiple time step algorithms Humphreys et al 1994 may be warranted The simulation of biological molecules at the molecular level has gen erated much excitement and these approaches have become an increasingly important partner with experimental studies of these complex systems 4 BIO INSPIRED MATERIALS Bio inspired materials represent a special area of opportunity for developing new high performance engineering materials based on ideas inferred from Nature Tirrell et al 1994 For example the proteins derived from spider silk serve as the inspiration for high strength fibers Simmons et al 1996 the adhesives from barnacles suggest how to produce glues that cure and function underwater and the complex protein inorganic interactions in mollusk shells supply ideas for producing ceramics that are less brittle than current ones It is likely that ultimate bio inspired materials will be chimeric that is they will be produced as a hybrid between biological and synthetic components Conse duently these materials represent a special class of the protein folding problem and of polymer physics II ORGANISMAL BIOLOGY The central organizing theme for Cells and Cell Systems is how behavior and function at one level of organization emerges from the structure and interactions of components at lower levels In the set of topics described in this section the lower level of orga nization is subcellular or cellular Though some of the subcellular components that play a role in these models are molecular the focus is not on the structure of those molecules but on the part that they play in cellular and multicellular function 1 CELL SIGNALING Control of cellular processes mediated by interactions of signaling molecules and their cell surface re ceptors is a central and unifying theme in current experimental 9 cell biology Within the past five years techniques of molecu lar biology have revealed many of the kinases phosphatases and other molecules involved in signal transduction pathways as well as molecular sub domains and sequence motifs that determine distinct functions New techniques for measuring phosphoryla tion calcium fluxes and other early biochemical responses to receptor interactions are being applied to study many cell sig naling systems eg chemotactic bacteria neurons and lympho cytes Genetically engineered experimental systems consisting of homogeneous cell lines transfected with homogeneous popu lations of wild type and mutant receptors and effector molecules have facilitated acquisition of much of the new information about the intracellular molecules that mediate signal transduction Im proved measurement and experimental design make mathemati cal modeling an increasingly feasible tool for testing ideas about the interactions of these molecules 2 MECHANICS AND EMBRYOLOGY Recent advances in instrumentation have made it possible to measure motions and mechanical forces at the molecular scale Svoboda and Block 1994 Concomitant with these new me chanical measurements are crystallographic and x ray diffraction techniques that have revealed the atomic structure and molecular geometry of mechanochemical enzymes to angstrom resolutions Rayment and Holden 1994 Together these techniques have begun to supply data that has revived interest in cellular mechan ics and reinvigorated the view of enzymes as mechanochemical devices It is now possible to make realistic models of molec ular mechanochemical processes that can be related directly to experimentally observable and controllable parameters Peskin and Oster 1995 These advances in experimental technology 10 have initiated a renaissance in theoretical efforts to readdress the central question how do protein machines work More precisely how is chemical energy transduced into directed me chanical forces that drive so many cellular events Embryology has also moved beyond descriptive observation to encompass genetic control of development and localization of protein effectors The stress and strain measurements that are now possible at the cellular scale promise to unite the genet ics biochemistry and biomechanics of development Oliver et al 1995 By characterizing the mechanical properties of embryonic cells and tissues mathematical models can be used to discrim inate between various possible mechanisms for driving morpho genesis Davidson et al 1995 3 BIOFLUID DYNAMICS Because of the ongoing revolution in computer technology we can now solve fluid dynamics problems in the three spatial di mensions and time Ellington and Pedley 1995 This opens up biological opportunities on many different scales of size On the organ scale for example one can now perform fluid dynamics simulations of the embryonic and fetal heart at different stages of development Such models will help to elucidate the role of fluid forces in shaping the developing heart The swimming mechanics of microorganisms are also accessible to computer simulation A particularly challenging problem in this field con cerns the intense hydrodynamic interaction among the different flagella of the same bacterium When the flagella are spinning so that their helical waves propagate away from the cell body 11 they wrap around each other to form a kind of superflagellum that propels the bacterium steadily along when their motors are reversed and the flagella spin the other way the superflagellum unravels and the bacterium tumbles in place Because of the difficulty of measuring microscopic fluid flows hydrodynamics within cells is a much neglected aspect of cellular and intracellu lar biomechanics Indeed computation provides our only window onto this important aspect of cellular physiology The incom pressibility and viscosity of water have the effect of coupling motions along different axes and between objects quite distant from one another biomolecular processes are also modulated by the necessity of moving water out of the way A new feature in this realm of micro and nano hydrodynamics is the importance of Brownian motion and the related significance of osmotic me chanics including sol gel transformations for controlling fluid motions 3 IMMUNOLOGY AND VIROLOGY During the last two years mathematical modeling has had a major impact on research in immunology and virology Serious collabo rations between theorists and experiment provided breakthroughs by viewing experiments in which AIDS patients were given po tent anti retroviral drugs as perturbations of a dynamical system Mathematical modeling combined with analysis of data obtained during drug clinical trials established for the first time that HIV is rapidly cleared from the body and that approximately 10 billion virus particles are produced daily Ho et al 1995 This work had tremendous impact on the AIDS community and has for the first time given them a quantitative picture of the disease pro cess The impact of this type of analysis has extended beyond AIDS and opportunities exist for developing realistic and useful 12 models of many viral diseases Challenges remain in studying drug therapy as a nonlinear control problem and the issue of how rapidly viruses mutate and become drug resistant under dif ferent therapeutic regimes needs to be considered Such issues also apply to the development of antibiotic resistance in bacterial disease 4 NEUROSCIENCES The fundamental challenge in neuroscience is to understand how behavior emerges from properties of neurons and networks of neurons Advances in experimental methodologies are provid ing detailed information on ionic channels their distribution over the dendritic and axonal membranes of cells their regulation by modulatory agents and the kinetics of synaptic interactions The development of fast computing sophisticated simulation tools and improved numerical algorithms has enabled the de velopment of detailed biophysically based computational models that reproduce the complex dynamic firing properties of neurons and networks Such computations provide a two fold opportu nity for advancing our knowledge 1 they both explain and drive new experiments 2 they provide the basis for new math ematical theories that enable one to obtain reduced models that retain the quantitative essence of the detailed models These reduced models which allow the bridging of multiple spatial and temporal scales are the building blocks for higher level models Modeling tools and mathematical analysis allow us to address the central question What are the cellular bases for neural compu tations and tasks such as sensory processing motor behavior and cognition Koch and Segev 1989 Bower 1992 More specifi cally how do intrinsic properties of neurons combine in networks with synaptic properties connectivity and the cable properties of dendrites to produce our interaction with the world Neu ral modulators affect both the intrinsic currents and the synap tic interactions between neurons Harris Warrick et al 1992 The effects of these changes at the network level are difficult to work out even for small networks The largest challenge in this area is to understand how systems with enormous numbers of degrees of freedom and large numbers of different modulators combine to produce flexible but stable behavior The geometry and electrical cable properties of the branching dendrites of neu rons also affect network activity Stuart and Sakmann 1994 Mathematical analysis is needed to interpret the results of mas sive computations and to incorporate the insights into network models III ECOLOGY AND EVOLUTIONARY BIOLOGY Evolution is the central organizing theme in biology eg Rough garden1979 and its manifestation in the relationships among types of organisms spans levels of organization and reaches out from biology to earth and social sciences Thus the core prob lems in ecology and evolution run the gamut from those that address fundamental biological issues to those that address the role of science in human affairs Fundamental challenges facing ecologists and evolutionary biologists relate to the threats of the loss of biological diversity global change and the search for a sustainable future as well as to the continued search for an un derstanding of the biological world and how it came to assume its present form To what extent is the organization of the biological world the predictable and unique playing out of the fundamental 13 rules governing its evolution and to what extent has it been con strained by historical accident How are the interactions among species ranging from the tight interdependence of host and par asite to the more diffuse connections among plant species in a forest manifested in their coevolutionary patterns and life his tory evolution What are the evolutionary relationships among closely related species in terms of their shared phylogenetic his tories How do human influences such as the use of antibiotics and pesticides exploitation of fisheries and land and accelerated patterns of global change influence the evolutionary dynamics of species and patterns of invasion To what extent can an evo lutionary perspective help us to prepare for the future in terms of understanding what species might be best suited to new en vironments The latter is important both in terms of natural patterns of change and deliberate manipulations through breed ing and species introductions 1 Conservation biology the preservation of biodiversity What factors maintain biodiversity How can new approaches to phylogenetic analyses in clarifying the evolutionary relationships within and among species help us to understand how we should measure biodiversity How are ecosystems organized into func tional groups ecologically and evolutionarily and how does that organization translate into the maintenance of critical ecosys tem processes such as productivity and biogeochemical cycles as well as climate mediation sequestering of toxicants and other issues of importance to human life on earth 14 2 Global change What are the connections between the physical and biological parts of the global biosphere and the multiple scales of space time and organizational complexity on which critical processes are played out Bolker et al 1995 In particular how are in dividual plants influenced by changes in atmospheric patterns and more difficult how do those effects on individual plants feed back to influence regional and global patterns of climate and biological diversity How do effects on phytoplankton and zooplankton relate to each other and to the broader patterns that may be observed 15 3 Emerging disease How do patterns of population growth and resource use as well as the profligate use of antibiotics contribute to the emergence and reemergence of deadly new diseases many of them antibiotic resistant Ewald 1995 Are there approaches to management of the diversity of those diseases guided by both an evolutionary and an ecosystem perspective that can reduce the threat and provide new strategies for mitigation 16 4 Resource management The history of the management of our sources of food and fiber is not one of unmitigated successes and many of these cru cial resources are threatened to a level that they will be unable to support the needs of humanity in the coming decades The prospect of large scale alterations of the earth s physical and biological systems creates a potential conflict between human needs desires and capabilities Walters and Parma 1996 Wal ters and Maguire 1996 This situation is further complicated by the limitations of our understanding and ability to control com plex biological systems We must develop methods for decision making and management that are appropriate for an uncertain future Hilborn et al 1995 17 A POPULATION GENETICS With the rapid accumulation of sequence data for entire genomes we are now poised to analyze the set of genes their order and or ganization codon usage etc across taxa Griffiths and Tavare 1996 and how and perhaps why this has evolved over time Thorne et al 1992 This requires an increased ability to model how information is represented and acted upon in biological sys tems Griffiths and Tavare 1996 based on tools from such fields as discrete mathematics combinatorics and formal languages Novel perhaps ad hoc formulations are needed to form the math ematical basis of genomic analyses because classical quantitative formulations of notions such as information similarity and clas sification all inextricably related to biology are inadequate Correspondingly methods for organizing vast sequence data into 18 data structures and databases suited for the most efficient data storage and access are needed along with improved algorithms for sequence analysis and the identification of homologies among sequences B CONSERVATION BIOLOGY Virtually all important questions in conservation biology require making predictions so theory and mathematical methods have played and will continue to play a central role Although many of the underlying scientific issues have been defined during the past decade many questions remain to be resolved What species would be lost in the wake of an invasion and what are the ef fects on ecosystem function For example what are the conse quences of the replacement of native fish species by introduced species Substantial progress is likely and needed in the near future in understanding the dynamics of invading exotic species determining more carefully the role genetics plays in the dynam ics of rare or endangered species and in the ecological dynamics of threatened species 19 C MANAGEMENT OF NATURAL SYSTEMS In recent years there has been an abrupt shift in management philosophy Hiborn et al 1995 The old goal of managing indi vidual species in order to reach and maintain optimal conditions has been replaced by a new goal of maintaining ecosystem func tion and adapting to new conditions or changes in the system This shift reflects a more mature attitude towards nature that recognizes the limitations of our knowledge and capabilities the importance of interactions between species and an appreciation of the dangers of a command and control mode of operation This new approach to management makes it possible to ap ply elements of the scientific method in a new and significant context we may design experimental management schemes to 20 provide information that is required to improve the management process and adapt to changes even unforeseen changes This new approach challenges our mathematical and statistical skills Successful adaptation reduires effective and timely organization of data through estimation of parameters that affect system dynamics including the dynamics of our learning That infor mation then must be translated into an assessment of the likely consequences of management strategies and actions D GLOBAL CHANGE AND BIODIVERSITY Climate change and associated changes in greenhouse gases have made imperative the examination of the potential impacts on natural systems and associated feedbacks Advances in compu tational capabilities have made possible the construction of de tailed individual based models that take account of the responses of individual trees to changes in environmental conditions and their mutual effects Yet such models are tremendously data hungry and have great potential for error propagation To make their predictions robust and to allow those predictions to be interfaced with the much broader scale predictions of climate models and the masses of broad scale information that are be coming available from remote sensing we must find ways to reduce dimensionality and simplify those overly detailed mod els Similar comments apply to models of other systems such 21 as the aggregation of social organisms from cellular slime molds to marine and terrestrial invertebrates and vertebrates Meth ods such as moment closure and hydrodynamic limits borrowed from other disciplines are proving remarkably promising espe cially when coupled with experimental approaches Levin and Pacala 1996 This represents one of the most challenging and important is sues in ecosystem science At the same time masses of data are becoming available from global observation systems and critical experiments are providing understanding of the linkages between ecosystem structure and function and in particular the role of biodiversity in maintaining system processes The next 5 10 years hold remarkable potential for integrated theoretical em pirical and computational approaches to elucidate profound and important issues Field 1992 Bolker 1995 E THE DYNAMICS OF INFECTIOUS DISEASES The subject of infectious disease dynamics has been one of the oldest and most successful in mathematical biology for a century and has seen powerful advances in recent years in mathemati cal theory and in the application of that theory to management strategies see for example Anderson and May 1991 Much of the literature has assumed homogeneous mixing so that every individual is equally likely to infect every other individual but such models are inadequate to describe the central qualitative features of most diseases especially those that are sexually trans mitted or for which spatial or socioeconomic structure localizes interactions The classical work of Hethcote and Yorke 1984 on core group dynamics highlighted the importance of such ef fects and formed the basis upon which much recent work rests Such work involving spatial structure frequency and density dependence and behavioral factors have not only forced us to revise old paradigms but have reenergized the interplay among nonlinear dynamics ecology and epidemiology Professor Donald Knuth on Bioinformatics Professor Donald Knuth of the Stanford Computer Science Department was once interviewed and had some comments on issues related to Bioinformatics CLB If you were a soon to graduate college senior or PhD and you didn t have any baggage what kind of research would you want to do Or would you even choose research again Knuth Ithink the most exciting computer research now is partly in robotics and partly in applications to biochemistry Robotics for example that s terrific Making devices that actually move around and communicate with each other Stanford has a big robotics lab now and our plan is for a new building that will have a hundred robots walking the corridors to stimulate the students It ll be two or three years until we move in to the building Just seeing robots there you ll think of neat projects These projects also 22 suggest a lot of good mathematical and theoretical questions And high level graphical tools there s a tremendous amount of great stuff in that area too Yeah I d love to do that only one life you know but CLB Why do you mention biochemistry Knuth There s millions and millions of unsolved problems Biology is so digital and incredibly complicated but incredibly useful The trouble with biology is that if you have to work as a biologist it s boring Your experiments take you three years and then one night the electricity goes off and all the things die You start over In computers we can create our own worlds Biologists deserve a lot of credit for being able to slug it through It is hard for me to say confidently that after fifty more years of explosive growth of computer science there will still be a lot of fascinating unsolved problems at peoples fingertips that it won t be pretty much working on refinements of well explored things Maybe all of the simple stuff and the really great stuff has been discovered It may not be true but I can t predict an unending growth I can t be as confident about computer science as I can about biology Biology easily has 500 years of exciting problems to work on it s at that level Graduate study in Mathematical Biology University of British Columbia http www math ubc capeoplefacultykeshetmathbiomathbio html RESEARCH AREAS INCLUDE Ecology Animal Behaviour Evo lutionary Genetics Cellular and Molecular Biology Excitable Cells Neurobiology Morphogenesis Pattern Formation Self Organization Nonlinear Dynamics Chaos FACULTY Leah Keshet Cytoskeleton and actin dynamics swarm ing and aggregation behaviour in animal societies Yue Xian Li Calcium dynamics signal transduction in cells biophysics neu roscience Michael Doebei mathematical models of evolution Fred Brauer Mathematical epidemiology and ecology 23 U n iversity of Washington http www amath washington edu RESEARCH AREAS Our interest lies in understanding the spa tial and temporal patterns that arise in dynamic biological sys tems Our mathematical activities range from reaction diffusion equations to nonlinear and chaotic dynamics to optimization We employ a variety of tools and models to study problems that arise in development epidemiology ecology resource manage ment and biomechanics FACULTY Mark Kot J D Murray H Qian 24 Princeton University httpwwweebprincetoneduindexhtml RESEARCH AREAS Ecology evolution and behavior FACU LTY Simon A Levin Dynamics of populations and communities spa tial heterogeneity and problems of scale evolutionary ecology theoretical and mathematical ecology biodiversity and ecosys tem processes Stephen W Pacala Population biology and community ecology of plants theoretical and mathematical ecology global interac tions among the biosphere atmosphere and hydrosphere 25 New York U niversity http www cns nyu edu RESEARCH AREAS The research interests of the faculty span a broad range of topics in neural science and utilize techniques ranging from molecular and cellular analyses to fully integrated systems computational and cognitive studies FACULTY John Rinzel Biophysical mechanisms and theoretical foundations of neural computations Charles S Peskin Mathe matical biology Michael J Shelley Professor of Mathematics and Neural Science Michael J Shelley Applied Mathematics Modeling and Large Scale Computation Vision and Computa tional Neuroscience 26 University of Pittsburgh http www math pitt edu G Bard Ermentrout Dynamical systems nonlinear differential and integral equations mathematical biology neuroscience Car son C Chow Applied mathematics biological mathematics non linear dynamics Jonathan Rubin W amp M graduate Dynamical systems computational neuroscience 27 Boston University httpmathbuedu Univ of Tennessee httpwwwmathutkedu Duke U niversity httpwwwmathdukeeducmclm NC State University httpwwwstatncsueduprogramsbmahtml University of Utah httpwwwmathutahedugradmathbiolindexhtml 29 Chapter 2 Diffusion equations 21 Separation of variables intervals Diffusion equation is a linear partial differential equation7 since the functions related to u in the equations ut and Au are both linear Recall for the linear ordinary differential equations 2 kP and DE kP 21 13 it is well known that 5 and Ac0s kth Bsin kDz are the solutions of the respective equations Indeed in the second case7 the exponential function skimD is a complex solution Thus we can make a lucky guess for the solution of the diffusion equation 8P 82F a Tm 2392 that the solution is an exponential function Pt z swim 23 for some constants 11 which could be complex numbers By substitution7 we nd that the function in 23 is a solution of 22 if a D1927 and depending on the values of b7 we obtain three families of solutions Ptz eDbZtebm b gt 0 24 Ptz eDV e bw b gt 0 25 Ptz 513172517 b gt 0 26 The rst two families of solutions are not reasonable solutions7 since when x a 00 or foo Pt7 z A 007 which implies an unlimited growth at z 00 or foo The third family of solutions are complex7 but their real and imaginary parts are both solutions of 22 too Pt x 5 13th cosbm and Pt x 5 13th sinbx 27 19 20 CHAPTER 2 DIFFUSION EQUATIONS These solutions are little more reasonable as they are bounded as x a 00 but still they do not satisfy natural boundary conditions P Pm a 0 as x a 00 Nevertheless solutions with above forms are solutions of diffusion equations and we notice that they are in form of Ptz UtVz This motivates us to consider the solutions which is separable on the two variables 25 and z utz UtVz 28 With the form of u in 28 the diffusion equation 22 becomes U t VW U WDU Wam Dk Ut 7 We 7 where k is a constant independent of both if and z Thus U satis es an equation U t DkUt 25gt 0 29 and Ut 015W 210 and V satis es V z 211 The solutions of 211 are Cge m with any complex number k Thus we obtain solution of 22 CeDkte m 212 just like the three families of solutions we obtain by making lucky guess To obtain more speci c solutions we must take the boundary and initial conditions into considerations We start with the one dimensional diffusion equation on an interval with length L 8P 7 82F 3t 7 8x2 7 Here the individuals of the the species inhabits a one dimensional patch 0 L and the diffusion coef cient of the population is D gt 0 To get a de nite solution of the problem we need to have an initial population distribution 35 e 0 L 213 P mean no and a boundary condition We assume that the exterior is hostile so homogeneous Dirichlet bound ary condition is satis ed Pt0 Pt L 0 215 We look for solution with form Pt z UtVz Then we obtain 29 and 211 and 215 implies that V0 VL 0 216 211 and 216 become a boundary value problem for an ordinary differential equation which is a special case of Sturm Lz39oum39lle problems We shall discuss this kind of problems in a separate section next 22 EIGENVALUES AND EIGENFUNCTIONS 21 22 Eigenvalues and eigenfunctions of boundary value problems We consider the Dirichlet boundary value problem 2 kg m e M ylt0gt yltLgt 0 217 We can see that 0 is always a solution of the equation So we shall look for solutions other than the trivial zero solution From elementary ODE knowledge we know that the general solution of y kg is one of the following three forms clam mm if k gt 0 218 01 022 if k 0 219 01 cos7kt 02 sin7kt if k lt O 220 One important characteristic of a boundary value problem is that it may not have a solution while the initial value problem always has a solution and the solution is unique if all terms in the equation are differentiable For 217 there is no solution when k gt 0 Indeed if there is a solution it must be yt cle t cge t for some constants 01 and 02 and we can solve the constants by 0 y0 cl 1 2 0 yL cle L Cge L But this system of equations is unsolvable if one gets 01 702 from the rst equation and plugs it into the second one one has e L LCz O and one must have 62 0 and 01 0 Thus 217 has only zero solution when k gt 0 It is also easy to check that 217 has only zero solution when k 0 So the only meaningful case is when k lt O 75 In this case 217 may still have only the zero solution but for some number k lt 0 it can have a solution yt which is not the zero solution Note that from the linear principle if yt is a solution of 217 so is 0 yt for any constant c and if y1t and y2t are solutions of 217 so is y1t y2t So for k which 217 has non zero solutions we call these numbers eigenvalues and the solutions of 217 eigenfunctions In fact historically the solutions of 217 were rst called eigenvalues and eigenfunctions Now let s search for the eigenvalues and eigenfunctions of 217 the eigenfunction must be yt 01 coshfist 62 sin7kt for some constants 01 and 02 by the boundary conditions we ave 0 y0 01 0 yL 01 cos7kL 6281H7kL so 01 0 and if 02 is not zero sinx7kL 0 thus xikL mm the eigenvalue k must be m27r2 k 77 221 2 t and the eigenfunctions associated with k imL are csin Therefore 217 has a mz 392 7 L2 sequence of eigenvalues km 7 and the set of eigenfunctions are the sine functions which 22 CHAPTER 2 DIFFUSION EQUATIONS can t in77 the interval 0 L When the interval is 0 7r so L 7139 the expression of eigen pairs is much simpler km 7m2 fmt sinmt 222 Back to 213 and 215 we nd a sequence of solutions 2 2 umtm exp lt7DmL t sin m is a positive integer 223 Because of the linear principle the general solution of 213 is 00 mz 392 7mm 142532 gcmemp 71372 s1n 25gt 0 m e 0 L 224 Example 21 y kg 0 220 we yltwgt 0 225 This problem is same as previous example except the boundary conditions So we still consider the cases of k lt O k 0 and k gt 0 When k lt 0 yt cle m Cge t from the boundary conditions we get y0 y 0 01 02 77 k01 x7 kcg 0 and y7r 1equot 1W 02 0 So 01 02 satis es the equations a wma i i i 1 7 7k 1 x 7k i The equation has non zero solution only if det lt sixj EH O The determinant equals 1 to 17 pe 7r 7 1 pe pquot if we denote p x7k thus det 0 is equivalent to 52pquot The 7 p i 2 r 1 p i 1 functions f1 p e p and f2p 1 have no common p01nts when p gt 0 in fact they only 7 p 1 i i i i p is not possible for p gt O and there is no eigenvalue such that intersect at p 0 thus 52 quot k lt 0 When k O yt 01 0225 from the boundary conditions we get y0 y 0 01 02 0 and y7r 01 027r 0 since the determinant of the matrix lt1 W is not zero then only solution for 01 02 is 00 So 0 is not an eigenvalue either For k gt O yt 01 cosh t 62 sinEt The boundary conditions y0 y 0 01 1 02 O y7r 01 coshEr 02 sinE7r O The matrix is MiEr sinkgwgt and the determinant of the matrix is 0 if sinE7r 7 xEcosE7r O or tan7rp p for p The graphs of f1p tan7rp and f2p p intersect at a sequence of points pm m 1 23 and 05 lt p1 lt 15 15 lt p2 lt 25 m 7 05 lt pm lt m 05 km p3 are the eigenvalues of the problem 225 and the eigenfunctions are ymt McosQEt 7 sinEt 23 FOURIER SERIES 23 23 Fourier series To completely solve 213 215 we need to determine the constants cm in 224 by setting t 0 in 224 00 u w 0 sin m 26 OL 227 olt gt m L lt gt lt gt The expression in the right hand side of 224 is the Fourier series of u0z on 0 L Fourier series is widely used in mathematical physics and signal process All sine and cosine functions are considered to be regular smooth waves and Fourier series convert any noise or signal into a sum of more regular waves and that is a big help for the signal processing Here we collect a few facts of Fourier series Theorem 22 Let and be piecewise continuous functions on the interval 0L Then can be expanded in either a pure sine series 102 bmsin 228 2 L where bm sin dz m 1 2 or a pure cosine series 0 102 22am cos 229 m n w 2 L where am ZJ fz cos lt L gtdxm012 The Fourier series is similar to the Taylor series n 1 2 f 090 230 m0 39 71 but instead of expanding the function in polynomials Fourier series expands a function in sine or cosine functions But Taylor series converges to the original function only in a small neighborhood of z 0 though some time it is large and the neighborhood depends on the function f Fourier series converges to the original functions for all functions and on the whole interval 0 L But it may not converge to fz at z 0 and z L Theorem 22 shows all sound wave can be decomposed into several regular waves plus a small noise The approximation of Fourier series is sine or cosine Fourier polynomials N N fzzmz1bmsinltmmgt fx amCOSlt m n w a 231 24 CHAPTER 2 DIFFUSION EQUATIONS Now we can complete the solution formula of 213 215 00 mz 392 7mm utz Z cmezp lt7D L2 25 s1nlt L gt quotF1 L 232 2 7mm gt where cm LO u0z s1n lt L gtdmm 71 25gt 07 z E 07 L Example 23 utummy tgt07 07r7 Mt7 0 Mt7 7r 07 233 140 z sin m From the solution formula 2327 we get 2 2 W utz Z cme m t sinmz7 and cm sinzsinmzdm 234 7quot 0 m1 7r 7r However sinzsinmzdz 0 for any m 2 27 and sin2 mdz 7r2 Therefore the solution is 0 simply 74257 z 5 sin m 235 Example 24 For the homogeneous Neumann boundary value problem ut Dum7 tgt 07 6 ELL7 72120 M057 L 07 236 M07 95 140907 the solution is 2 0 2 74257 z 630 Z cmezp lt7DmL t cos m1 77mm 2 L where cm Zj u0z cos lt L 237 gtdz mZO tgt0 x 0L The derivation of this formula is left as exercise For these examples7 the asymptotic behavior of the solution can be read from the formulas For homogeneous Dirichlet problem 213 2157 the asymptotic limit is i i i 7r2t i 7m lim utz O7 and the asymptotic pro le is utz N clamp 7 2 sin 238 taco L L On the other hand7 for homogeneous Neumann problem 2367 the asymptotic limit is 2t 74257 z 7 and the asymptotic pro le is 74257 z N 630 clamp cos 239 24 RELAXATION TO EQUILIBRIUM SOLUTIONS 25 L The constant co above is 2L u0zdm thus the limit is 0 L A mam 240 the average value of the initial value u0m The limits of both cases are also solutions of diffusion equation and in fact they are solutions which are independent of time 25 Such solutions are equilibrium solutions For a general reaction diffusion equation 8 8 dAufmu tgt0 m E 9 Mom 140907 m e a 241 Bu0 tgt0m 8 2 where Bu is the boundary condition the equilibrium solutions satisfy dAu 7 0 m E Q 242 Bu 0 m E 89 A solution u1z of 242 is a solution of 241 in the sense that the function U1t z u1z is a solution of 241 with u0 z It is not just a coincidence that for two examples we have considered so far the limit is always an equilibrium solution We will show in next section that the limit of a diffusion equation is usually an equilibrium solution Equilibrium solutions of diffusion equation are easier to nd and that will help us understand the limit of solution when the solution itself is not so easy to nd The equilibrium equation for 213 215 is u 0 u0 uL 0 243 The general solution for u 0 is ax b and from the boundary conditions we obtain 0 is the only equilibrium solution The equilibrium equation for 236 is u 0 uO uL 0 244 For this problem b for any I E R is an equilibrium solution However the limit in 236 is the average value of the initial value u0z only Why does the solution of diffusion equation selects this particular equilibrium solution but not the others We will answer this question in next section 24 Relaxation to equilibrium solutions We have obtained analytic solutions of 1 d diffusion equation in previous sections The solutions are in a form of in nite series However in general even such formulas cannot be obtained when the spatial dimension is higher and the domain is arbitrary Later we will formulas of solutions of diffusion equations for some special cases Here we use some integral formulas to show some qualitative properties of solutions of general diffusion equation 26 CHAPTER 2 DIFFUSION EQUATIONS As a calculus warmup we consider the diffusion equation with no ux boundary condition jag 7DAu tgt0z 2 74071 was x e o 245 Vutm nx 0 t gt 0 m E 89 From the biological modelling point of view the species only randomly moves in Q the growth rate is zero and the boundary is sealed Thus the total population should not change since no one moves in or out and no one is born or dead lndeed that is the basic character of the equation and it also explains why the limit of the solution of 236 is the average value of the initial value function Proposition 25 Suppose that utm is the solution of 245 then utzdz u0mdm 246 9 9 Proof Integrating the left hand side of the equation in 245 over 9 we obtain 9 tmdm ltQutmdmgt 247 and integrating the right hand side we have D Autzdm D was x nmds o 248 0 39 from the divergent theorem and the no ux boundary condition Therefore from 247 and 248 we nd that d E utxdmgt o 249 which implies 246 III Next we show that the solution ut z always tends to a limit constant c and because of 246 that constant must be 1 uo m dx 250 m 9 lt gt lt gt A rigorous proof of this fact needs many more advanced mathematical knowledge so we only sketch the ideas here We multiply the equation in 245 by the function ut m 1314 U5 7 DuAu 251 Similar to the proof of Proposition 25 we integrate the left hand side of 251 ea 7 a 12 7d 12 71d 2 Quatdmi at lt2ugtdmidt Qlt2ugtdmi2dt lutml dx 252 24 RELAXATION TO EQUILIBRIUM SOLUTIONS 27 and the right hand side D uAudx D uVu nds 7 D Vu Vudx 253 D 39 9 Here we use Green s identity a consequence of divergence theorem uAudx uVu nds 7 Vu Vudx 254 D 39 D From the boundary condition fanuVu nds 0 thus from 252 and 253 we nd 1 d ut m2dm 7D qut ml2dz 255 26h Q 9 Let Et ut m2dz Then 255 shows that E t lt 0 for all t gt 0 and obviously Et gt 0 Q for all t gt 0 Thus limtH00 Et 2 0 exists and as t a 00 E t a 0 so i 2 7 qutml dm 7 0 256 256 implies limtH00 qut 0 for any x E 9 thus the limit of utz when t a 00 must be a function with zero gradient thus the limit function must be a constant 0 Since elm tlim utzdm u0zdm then c Notice that we 9 4 00 Q Q implicitly assume utz has a limit when t a 00 While such a limit indeed exists the proof of its existence is beyond the scope of this notes that is why this is only a sketch of proof Nevertheless we put our ndings into another proposition Proposition 26 Suppose that utm is the solution of 245 then utm tends to a constant xdx function W uo 9 Similar result can also be shown for Dirichlet problem in which the limit function is u O From this property of diffusion equation we can get a conclusion that the diffusion equation will make the initial value atter and bring it eventually to an equilibrium state without spatial pattern This can also be indicated from the differential equation itself Let the average of the initial value be c For the points where ut z gt c the function is likely to be convex down with Au lt 0 thus the function value will decrease over these points and for the points where ut z lt c the opposite occurs This is called smothering effect of diffusion In the following we use a Maple program to demonstrate the smothering effect and the asymptotic tendency to the equilibrium We consider ut 4um t gt 0 m 6 04 ut0 ut 4 0 257 u0z 28 CHAPTER 2 DIFFUSION EQUATIONS 100 x 6 13 From 232 we get the 0 x E 01U34l where x is given by the formula x formula of the solution of 257 7 00 m27r2 m7rx utx 7 gcmexp lt7 4 t sin 4 258 1 where cm 5 xsin dxm 21 25gt 0 x 6 04 0 The coef cients cm in the Fourier series can easily be calculated but we will let Maple to do that Here is the Maple program diffusionmws and some of the results gt restart gt withplots n n gt fx gt100Heavisidexl100Heavisidex3 gt plotfx x04 gt cn gtO5intfxsinnPix4x04 Cn cos14 n Pi cos34 n Pi 2000 n Pi gtuapprox1xt gtsumcnsinnPix4exp4nquot2t n1 30 30 uapprox1 xt gt Z cn sin14 n Pi x exp 4 n2 t n1 gt tvals seq0 0720i i0 20 gt toplotsequapprox1xt ttvals gt plottoplotx0 4 gt densityplotuapprox1x t x0 4 t0 O 2 colorstyleHUE gt animateuapprox1xt x0 4 t0 O 1 frames50 This is an animation generated by approximating the solution by the partial sum with n 30 of the series solution The approximation is accurate since the remainder is an exponential decaying function with smaller exponent Note that the initial condition x can be written as the difference of two Heaviside functions which is de ned as 1 x Z 0 0 m lt 0 259 It is an interesting experiment to observe the changes of the spatial patterns of ut x when 25 gets forward Here are a few snap shots of the animation 24 RELAXATION TO EQUILIBRIUM SOLUTIONS 29 Figure 21 a fm b partial sum of Fourier series of fz 100 100 80 80 so so 40 40 2o 20 0 1 g 3 4 0 1 g 3 4 Figure 22 a u0002m b u0010z Because of the truncation of the series solution utz above indeed is the solution of 257 with initial data 2 cm sin dm cm cos 7 cos 260 4 m H Figure 21 is a comparison of f and From the last section we know that when t a 00 ut z a 0 exponentially fast However the road to distinction is worth a more careful lookingiit can be dismantled into the following stages 1 Forming a plateau From 25 0 Figure 21b tot 0002 Figure 22 a we can observe that the initial data has been smoothed by the diffusion equation and all the small wiggles in u0z are gone at t 0002 On the other hand although u0 z is not a constant for z 6 13 u0002z is almost a constant on an interval slightly smaller than 13 So in the initial stage of the evolution of 1 under diffusion a plateau or at core of height 100 forms with a sharp but smooth drop off to 0 near z 1 and z 3 CHAPTER 2 DIFFUSION EQUATIONS Figure 2 03 q Figure 24 a u0100m b u0600z Erosion of the plateau From 25 0002 Figure 22 a tot 001 Figure 22 b the at core keeps shrinking though the top of the core is still at the level near 100 At t 001 the erosion nally dissolves the whole at top and the maximum point of ut at z 2 looks more like a regular local maximum point The graph of u001 z is close to a bell shape In that time interval the in ection points of ut are still near z 1 and z 3 but at t 001 the interface between u 0 and u 100 is no longer as sharp as when t 0002 Bell getting round In the next stage smoothing effect brings the concave part of the graph down and the convex part of the graph up And at the same time the convex part is shrinking as population loss via boundary gets bigger At t 004 the graph becomes almost concave for all z and the bell shape becomes an arch see Figure 23 Collapse of the arch In the nal stage the arch is close to the leading term of the series clezp77r2t4 sin7rz4 and an exponential collapse is obvious After 25 1 the graph can hardly be seen see Figure 24 25 ONE DIMENSIONAL CHEMICAL MIXING PROBLEM 31 25 One dimensional chemical mixing problem Here we consider a simpli ed version of the chemical mixing problem which is discussed in Chapter 1 A tube with length L contains salt water The cross section yz direction of the tube is so small so we can assume that the concentration of the salt water is same for any point on a cross section Let the concentration of the salt water in the tube be ctm kgms t gt 0 and 0 lt z lt L We assume that the diffusion of the salt is onedimensional in the tube and the velocity of uid is ignorable thus 005 z satis es the diffusion equation 36 82 E 261 At the left hand side of the tube a salt water solution with constant concentration C0 kgm3 enters the tube at a rate of R0 77138 and on the right hand side of the tube the mixed solution is removed at the rate of R0 mSs We also assume that the area of the cross section of the tube is A Let A1 and A2 be the cross sections at z 0 and z L correspondingly From the Fick s law the total in ux at z 0 is J inds 713mg 0 1ds iDAcmt0 00120 262 A1 A1 and the total out ux at z L is J nds 713mg L 1ds iDAcmt 0 Ct LR0 263 A2 A1 thus we obtain the boundary conditions 8675 0 CoRo 8675 L 7 7675 LR0 7 d 7 264 8m DA an 8m DA To simplify the notations we de ne 7 CoRO 7 R0 B7 DAandE7DA 265 Thus the initial boundary value problem is 2 D tgt0 ze 0L 13 80t 0 8025 L 2 66 7B E t L 0 39 8 8 clt gt ctmcoz tgt0 m 0L Before attempting to solve the equation we could nd the equilibrium solutions which satisfy equation 82 138 0 m e 0L 13 860 7 B 80L 8m 8m 267 EcL 0 32 CHAPTER 2 DIFFUSION EQUATIONS From a simple calculation we nd the unique equilibrium solution 7 B 7 CORO 01m7EBL7m 700 which is a linear positive function with negative slope on 0L and 01L Co To solve the equation we de ne bt z Ct z 7 01z Then bt z satis es 2 g 8 1 tgt0 0L L 0 0 L L Ect L o bt 00m 7 ciz tgt 0 m E 0 L L 7 x 268 269 The reason we consider 269 instead of 266 is that the boundary conditions in 269 is homoge neous thus we can use the technique of separation of variables We assume that bt z UtVz Then Ut V W DVz Dk 270 Ut 5D and Vm satis es V 7 kV 0 m e 0L V 0 0 V L EVL 0 271 This is a half Neumann and half Robin boundary condition We use the method in Section 22 to solve 271 If k gt 0 then Vm 1equotEm cge m From the boundary conditions we get 77 7E detlt Ef V L E e L gt 0 272 and k must satisfy 27H 7 E 5 7 xE E But the function on the left is greater than 1 when k gt O and the one on the right is always less than 1 when k gt 0 So 273 has no solution thus there are no positive eigenvalues Similarly k 0 is not eigenvalue either 273 If k lt 0 Vm 01 cosx7km 2 sinx7km From the boundary conditions we get 2 0 and 01cosx7kL 7 x7ksinx7kL 0 Thus k satis es cot7 kL 77 17 274 From the graph of the functions f1 cotLm and f2z 2 there are exactly one intersection point mm of f1 and f2 in m717rLm7rL for any positive integer m gt 0 Then this km 7z is the m th eigenvalue of 271 and the corresponding eigenfunction is cos Wm Thus the solution of 269 is m717r M 00 bt z Z bme Dmfnt cosmmz where mm 6 L L satis es cotsz mm 275 m1 26 CRITICAL PATCH SIZE AND BlFURCATIONS 33 Hence the solution is the sum of the equilibrium solution and an exponential decaying Fourier series and the solution converges to the equilibrium solution exponentially fast It shows that in equilibrium state the exiting solution also has the concentration CO same as the incoming solution thus the amounts of salt entering and exiting the tube are balanced but the concentration of solution in the tube is higher than CO and the concentration at z 0 is the highest The example in this section can be thought as a model of a continuously polluted river The left end point z 0 is the source of the pollution where the pollutant enters the river at a constant rate and the right end point z L is where the river merges to a bigger river or the ocean Our mathematical result shows that the concentration of pollutant in a continuously polluted river is even higher than that in the source As a consequence we can see that at the con uence of two polluted rivers the pollution can be most serious since the concentration there is much higher than the ones of either branches We should be cautious that we ignore the effect of drifting here which could make a big difference for rivers 26 Critical patch size and bifurcations The simplest population growth model is the Malthus model dPdt aP for a positive constant a gt 0 In this section we consider Dirichlet boundary value problem of diffusive Malthus equation 2 8 13 Da PaP tgt0 m 0L at 8x2 ut 0 ut L 0 276 u0x u0x m E 0 L We have seen that population are at a loss from the hostile boundary but on the other hand the population will increase through linear reproduction So the loss and gain are in a competition Let s see what factor will determine the outcome of this battle Equation 276 can be solved in a similar way as 213 215 In fact the form of the solution only changes slightly adding a term 5 to re ect the exponential growth The solution of 276 is co 2 2 170590 Z cmezp at 7 DmL t sin 7 m1 L where cm u0zsin dmm 21 25gt 0 m e 0 L 0 277 However the fate of the solution may be different since at is added to the exponent now the exponent may be a positive one then ut x will have a part which is exponentially increasing In fact when the index m in the sum increases the exponent gets more negative thus whether ut z is exponentially increasing is determined by the rst exponent coef cient a 7 D39n39ZLZ We have D712 D i whenagt 7 orLgt 471 tlimou m foo 2 278 D7r D i whenalt70rLltHZ7n tlimoutz70 34 CHAPTER 2 DIFFUSION EQUATIONS If we assume that the growth rate a and D are intrinsic parameters determined by the nature of species and environment then the size of the habitat L will play the deciding role here The number L0 Dmr is called the critical patch size as the population cannot survive if the habitat is too small the outgoing ux wins over the growth but the population will not only survive but thrive to an exponential growth if the habitat is large enough growth outpaces the emigration We say that a bifurcation occurs at L L0 as the equilibrium solution u 0 changes stability type from stable for L lt L0 to unstable for L gt L0 We notice that at L L0 the equation has other equilibrium solutions sin7rmL the eigenfunction associated with rst eigenvalue k1 77r2L2 We will have more discussion on bifurcation problems when the growth rate is nonlinear From 278 we can also use a or D as bifurcation parameter instead of L sometime that is more convenient since we do not need to change the domain but we get equivalent results The concept of critical patch size is related to habitat fragmentation in ecological studies Habitat fragmentation1 is the breaking up of a continuous habitat ecosystem or land use type into smaller fragments which is considered to be one of several spatial processes in land transfor mation It is commonly used in relation to the fragmentation of forests Habitat fragmentation is mainly caused by human activities such as logging conversion of forests into agricultural areas and suburbanization but can also be caused by natural processes such as re From our mathematical result in this section if the the fragmentation of habitat limits the spatial movement of the plants and animals then the species can become extinct 27 Separation of variables rectangles Most of this chapter we have exclusively considered one dimensional spatial domain for the sim plicity of the mathematical analysis However any interesting application happens in 2 D or 3 D spaces In this section we discuss the separation of variables in a two dimensional rectangle We rst consider a linear diffusion reaction equation 8U 82 LL 82 LL a ltW8 y2gtu7 tgt07 79ER07 1X0757 Vun0 my 8R Mommy uo7y7 any 6 R 279 where D a b gt O and 8B is the boundary of rectangle R From the method of separation of variables we assume that ut my UtVmy 280 Then similar to onedimensional case we obtain U t Dk Ut 281 1From TEMS Terrestrial Ecosystem Monitoring Sites httpwwwfaoorggtostems 27 SEPARATION OF VARIABLES RECTANGLES 35 and 82 82 V V kV 25 0 R 0 0 1 WWW gt we mm gt 2827 VVn07 x734 68R For 2827 we need a further separation of variables Waxy WZy7 283 Then we obtain that W77 Z77 z y 2847 k7 WW Zy so both and Z yZy must be constants On the other hand7 the Neumann boundary condition implies that W 0 W a 0 and Z 0 Z b 0 Therefore W and Z satisfy W x k1Wm7 m E 07a7 WO Wa 07 285 Z y k2Zy7 y E 075 ZO Zb 07 286 and k k1 kg 287 Both 285 and 286 are familiar eigenvalue problems in one dimensional space which have been studied in Chapter 2 The eigenvalues and eigenfunctions of 285 are 774271392 n n39m km 777 27 W748 cos 7 n 071727 7 288 and the eigenvalues and eigenfunctions of 286 are m27r2 2 m n39y kzmi 7Zmycoslt m071727 289 Therefore the eigenvalues and eigenfunctions of 282 are 71271392 m27r2 kmm 7 2 7 b2 7 Vnmz7y cos cos 7 290 where n7m O7 17 27 Thus we nd that the solution of 279 is ut7z7y i cmmemkmw w cos cos 7 291 n0m0 where 0mm can be determined by the initial conditions The spatial patterns of the eigenfunctions for rectangle are much more complex than those of onedimensional Here we consider a square R 077139 x 077r7 then the rst ve distinctive eigenvalues are see the contour graph of eigenfunctions below k00 07 k01 k10 17 km 27 k20 02 47 k2 k1 57 292 36 CHAPTER 2 DIFFUSION EQUATIONS Figure 25 From left to right Figure 26 From left to right a V0 cos2y b V10 cos2z c V12 cosz cos2y d V11 cos2z cosy We notice that many eigenvalues km have multiplicity more than 1 for example km kmm 7m2 712 In such case the eigenspace is of higher dimensional so more possible spatial patterns are generated For example km k0gt1 71 and the eigenspace is spancosmcosy thus k1 cosz k2 cosy is a spatial pattern for any k1 and k2 In particular the patterns generated by cosz i cosy are symmetric with respect to one of the diagonal lines Figure 27 From left to right a cosz cos y b cosy 7 cos x The solution 291 is always an exponential growth if A gt 0 Indeed from the equation the no ux boundary condition prevents the emigration through the boundary and gt 0 implies a positive growth rate thus the total population in R has an exponential growth We can consider corresponding Dirichlet boundary problem 8 82 82 8 7 Wye QM tgt07ltz7ygteRlto7agtxlto7bgt7 Mammy 07 any 68R u0my u0zy my 6 R 293 28 STEADY STATES IN RADIALLY SYIVIMETRIC DOMAINS AND TRANSIT TIMES 37 Similar to the method above the solution is 00 71713913 m7r ut x y Z bmmeltkamAgtt s1n s1n 294 n1m1 where bmm can be determined by the initial conditions and kmm is de ned in 290 Similar to last section the persistenceextinction of the population can be determined when A gt 1 1 lim t u m oo D7r2 a2 1227 taco 7 7 2 95 1 i when W lt a Z g tlimou m 7 0 28 Steady states in radially symmetric domains and transit times Chapter 2 Exercises 1 Find the and f t of Neumann boundary problem y kg m E 0 L yO yL 0 296 2 Find the eigenvalues and eigenfunctions of periodic boundary problem 2 kg m 6 lt0 7r 210 m we W 297 3 Find the eigenvalues and eigenfunctions of Robin boundary problem I gt 0 y kg m 6 lt0 L we 12240 2L 42m 298 4 Find the eigenvalues and eigenfunctions of the problem yH2yky0 m6 07r y0 y7r 0 299 5 ln Example 21 use Maple to nd the numerical values of the rst three eigenvalues and plot the graphs of the corresponding eigenfunctions a Show that the solution of ut Du tgt 0 m E 0L umt 0 umt L 0 2100 M0790 14000 z 6 0L t cos L where cm u0mcos dm m 2 0 tgt 0 z E 0 L 0 co co mzn392 utm 3 Z cmezp lt7D L2 m1 2101 to H H H H to CHAPTER 2 DIFFUSION EQUATIONS a Find the solution of ut Alum tgt 0 m 6 01 um 0 Mt 1 0 2102 1407 90 1 05 1 m e 0 05 71 m e 051 b Modify the Maple program in Section 24 to simulate the solution of 2102 where fz is given by the formula fz Find the solutions of equilibrium equation u 0 m 6 01 u0 310 5 u1 7 511 7 2103 Find the solutions of equilibrium equation u 31 0 m 6 01 u0 5 u1 4 2104 Suppose that ut z is the solution of 2 7DAu tgt0 69 Mow was as e 9 2105 Vutm nz a utm 0 tgt 0 m E 89 where a gt 0 Show that the mean squared value 25 ut z2dm is strictly decreasing 9 unless ut z E 0 Suppose that ut z is the solution of the convective diffusion equation D gt 0 and V gt 0 814 8214 814 a Dy 7 Va m E 07r u0x u0x m E 07r 2106 u 814 020 7 t r 7 0 1 Show that the total population ut mdz is a constant if u0z sinz but not a constant if u0z cosz Suppose that in the chemical mixing problem in Section 25 the salt solution drifts to the right with velocity V Then Ct z satis es the convectivediffusion equation 36 82 36 EDW7VE 2107 with DV gt 0 Suppose that ctm satis es the same initial and boundary conditions in 266 2 00 H OJ H q H U H a STEADY STATES IN RADIALLY SYIVIMETRIC DOMAINS AND TRANSIT TIMES 39 a Find the equilibrium solution of 2107 b Show that the equilibrium solution is positive7 decreasing and concave down7 and explain the effect of the drifting to the geometric properties of the equilibrium solution c Compare the exiting ux 7D0mt L of the equilibrium solution with that of 267 Suppose that in the chemical mixing problem in Section 257 the salt is consumed in a chemical reaction at a rate k Then Ct7 z satis es a linear diffusion equation D82C 7k 815 3m2 3 2108 with Dk gt 0 Suppose that Ct7 satis es the same initial and boundary conditions in 266 a Find the equilibrium solution of 2108 You can use Maple to nd the general formula b Show that the equilibrium solution is positive7 decreasing and concave up7 and explain the effect of the reaction to the geometric properties of the equilibrium solution c Compare the exiting ux 7D0mt L of the equilibrium solution with that of 267 d Plot the graphs of the equilibrium solutions of 2677 2107 and 2108 with same parameters in Maple Determine the critical patch for a diffusive Malthus equation with Robin boundary conditions 2 2 1 DZ IaP7 tgt07 m 0L7 1 71020 We mm ebulttLgt 09gt u0x 7Lx7 m E 07 L Consider the doubly periodic boundary value problem 814 8214 8214 a D U7 tgt07 79ER07 1X0757 M0721 Why 14199 Mm y 2110 77x7 0 77x7 b7 Uy7 0 Uy7b7 u0my ugxy7 my 6 R Find the eigenvalues and eigenfunctions of the problem7 and the series representation of the solution Consider a convectivediffusion equation 814 7 8214 814 7 4 8 t 0 8t 3m2 825 7 quot gt 6 7Lt7 0 07 7Lt7 L 07 u0x 7Lx7 m E 07 L 07L7 2111 Find the solution of the equation with the following steps Chapter 3 Diffusion of point source and biological dispersal 31 Diffusion of a point source Suppose that a new species is introduced to an existing environment then naturally the individuals of the species will disperse from the release point That would similar to the diffusion when a drop of ink enter a jar of pure water Reaction diffusion equation can be used to describe an ideal growth and spatial diffusion phenomena The simplest way is to consider the dispersion in an in nite domainiwe rst consider the simplest one dimensional problem Suppose that Ptm is the population density function of this species 25 2 0 and z E foo 00 Suppose that at t O we import a group of M individuals to a point say z 0 Let u0z be the initial distribution of the population Then 00 imam M 31 700 However u0z 0 for any x 7 0 since before the dispersion the species does not exist in the environment yet On the other hand the population density at z 0 is given by the following limit fix u0xdx 00 311 T 0039 32 Thus u0z is an unusual function can be characterized by f 0 u0z 0 1 a f and imam M 33 0 otherw1se 700 This function is called a 6 function delta function in physics and we denote it by 60z a function which is in nity at a would be 6a Although this function is not a function in the normal sense we can check that the function M lt1gt t m ef 34 rth 37 38 CHAPTER 3 POINT SOURCE is a solution of it Du tgt 0 m E R u0x M60m 35 in the sense that I satis es the equation when t gt 0 and z E R and l ltIgtt 0 0 l ltIgtt0 36 Agghz z7tig17oo The function I is called the fundamental solution of onedimensional diffusion equation For each xed t gt O ltIgtt is a normal distribution function thus the graph of gtt is a bell shaped curve Moreover since the individuals in the population only move around in the environment there is no new reproduction the total population keeps constant 00 ltIgtt dx M 37 Biologically dispersion in a two dimensional space is more practical Similar to onedimensional case the equation Ut DUzm l Win7 75gt 07 7 39 6 R27 u07 7 39 600 7 3907 3398 has solution M 2 2 ltIgttmy 6333 39 Here 600m y is a function de ned as 00 if my 00 600907y 0 and R2 600907 yldmdy M 310 otherwise We also notice that ltIgtt z y is radially symmetric thus the diffusion from a point source in two dimensional space is radial which can be regarded as a consequence of the rotational symmetry of the Laplacian operator 32 Mathematical derivation of the fundamental solutions In this section we will solve the equation ut Du t gt 0 m E R u0z 60z 311 and derive the formula of the fundamental solution 34 One can see that the fundamental solution is not a function in form of UtVz thus the method of separation of variables won t work here But it also shows that there are other types of solutions of diffusion equation which cannot be obtained by separation of variables When deriving the diffusion equation in Chapter 1 via random walk an important assumption we make is that Az22At a D when the step sizes Am and Am approach zero So we can guess that that the scaling zZt can be a factor in the solution of diffusion equation Indeed if ut z is a 32 MATHEMATICAL DERIVATION OF THE FUNDAMENTAL SOLUTIONS 39 solution to the diffusion equation so is ua2t am a dilation of ut Thus we try a test function ut z vm2t and hope it may solve the diffusion equation Then 1 1 3 3 05 4g 11372 g 312 Substituting 312 into 311 and after some algebra we have 2D 1 y U 3 4D 3 0 313 By integrating we nd that My cly lZe y Dl 314 thus the general solution of 313 is 1134 cl y 27125714Ddz 2 315 0 for constant 01 62 E R So we now have another family of solutions of diffusion equation z2t utm 11m2t 31 zingTzLxmdz 2 316 0 However this solution is still not the fundamental solution We notice that if ut z is a solution of the diffusion so is its partial derivatives as the equation is linear We differentiate 316 with respect to m and we obtain 2 mm 2ei 317 Finally we choose 01 M47rD so that fR vt mdz M and we obtain the formula of the fundamental solution There are at least another three ways of deriving the fundamental solutions One is to use a test function ut z t vt m with unspeci ed Oz and B the second is to use Fourier transform an integral transform and the third is to use random walk approach in Chapter 1 and the central limit theorem in probability theory The solutions given in 316 are also very useful in some cases By making a change of variables 2 4Dw2 and assuming that ut z a 0 as x a ioo we can rewrite the solutions in 316 as 2 utm Cg lt1 7 4 57w dwgt 1317 erfmm 318 0 where erfy e wzdw is called the error function which is w1dely used in statistics and 7quot 0 applied mathematics The solution in 318 satis es 814 8214 D t 0 0 8t 8352 gt dim Ut 0 Cg lim utm 07 319 mace u0x 0 m E 0 The equation 319 is a diffusion on a half line with continuous source instead of a point source Here x 0 is the source where the concentration or density of u is kept as a constant 03 40 CHAPTER 3 POINT SOURCE 33 Chemical Pollution revisit In Section 25 we consider a chemical mixing diffusion problem in a tube with nite length L and a continuous source ux at one end Here we apply some results in two previous sections to consider some more chemical questions Example 31 Prof Shi s house is on a long boat canal One day a neighbor 200 meters away from his house has a small fuel spill 1 kg Suppose that the current in the canal is negligible and the fuel is only transported through diffusion Assume that fuel mixes rapidly across the width and stays on the surface of the canal and assume the canal is 10 meter wide and the diffusion constant D 0047712 s H How long does it take for the spilled fuel to reach Prof Shi s house to What is the concentration of the fuel at Prof Shi s house at that time 03 When does the fuel achieve the maximum concentration at Prof Shi s house What is the maximum concentration Let ft x be the concentration of the fuel Then 1 satis es the diffusion equation and the fuel spill can be thought as a point source release The long canal can be assumed to be in nitely long Thus 1 satis es ft Dfm tgt 0 m E R u0x 60z 320 The solution is readily given by the fundamental solution 302 1 ltIgttz 7 We 4Dt 321 However to give a meaningful answer to the question we have to de ne what means the fuel reaches the house In fact for any 25 gt 0 ltIgtt z gt 0 for all z E R so the fuel technically reaches the house immediately though only very tiny portion at a point far away from the source One way to de ne a minimum traceable level k and to calculate the time when the concentration at the given point reaches that level But a conventional way is to de ne the edge of the diffusing patch of the fuel We notice that ltIgtt z is a normal distribution function with variance 0 2Dt The patch is usually de ned as the spatial domain within i2a of the center Thus the edge is at z 22Dt and the time that the edge reaches the house can be calculated by 2V2Dt 200 322 and t 1068 1250003 z 32 hours The concentration at that time is 1 7 10x47r 004 125000 The time 32 hours is when the edge of the patch reaches the house but not the time the maximum concentration is at the house To calculate the maximum concentration we can use Maple to differentiate the solution with respect to t and solve the time of maximum concentration 25 125000 2x2Dt 5 2 54 x 105kgm2 0054gm2 323 34 DIFFUSIVE MALTHUS EQUATION AND BIOLOGICAL DISPERSION 41 5000003 z 128 hours The maximum concentration is about 0121gm2 We also notice that although when t a 00 the solution tends to zero but the rate of approaching zero is extremely slow For instant the time of ft 200 back to 0054gm2 is 63021153 x 1750 hours Thus even the peak has passed after 5 days of the spill there is still a noticeable trace of the fuel in the period of 73 days Now you can imagine how bad it would be if a whole fuel tank spills on a highway or a oil vessel spills in the ocean u arm n mm D uma n mm D mm D mm D zumuu 40mm amuuu anan was Figure 31 The graph of ft 200 34 Diffusive Malthus equation and biological dispersion The simplest population growth model is due to Thomas Malthus 1P 7 dt 7 where a is the constant growth rate per capita Thus if we take the spatial distribution into the consideration it is natural to consider the following diffusive Malthus equation in two dimensional space aP 324 8P 82F 82F P 325 a 8x28y2a7 where a is the constant growth rate per capita per unit area In this section we consider the equation in the whole space with a point source 8P 82F 82F 2 EiDW8 y2aP tgt0 my R 74079079 600957 326 The solution of 326 can easily deduced from the solution of 38 the 2 D fundamental solution In fact let ltIgtt m y be the solution of 38 then it is easy to verify that eatltIgtt z y is the solution of 326 Thus the solution of 326 can be expressed as w2y2 M Ptzy 47rDtgat 4Dt 7 327 42 CHAPTER 3 POINT SOURCE and consequently Pt z ydzdy Me 328 R2 This equation can be used to model the population explosion of modern metropolitan area in twentieth century Suppose that A is an isolated metropolitan area which is far away from any other metropolitan areas then we can assume that A and the surrounding areas is an unlimited area RZ Here we also assume that the migration of the population can be in any directioni there is no water bodies in the area no unlivable mountain desert etc At the time t 0 the beginning of the city a group of M immigrants came to this place and began the formation of the city by diffusion and growth After some time when a town has been formed near the center of the rst settlement 00 we can de ne the metro area and the rural area in the following way Let BR be the ball with radius R in R2 we de ne Rt as the number such that outside of BRm the total population is always M 239e Ptzydmdy M 329 RZ BRO we de ne BRm as the metro area and R2 7 BRm as the rural area We make such a de nition in assuming that the population in all rural area keeps a constant despite the growth but the population in metro area has an exponential growth and the area of the metro also grows over years With the calculation which leads to 328 we have R2t Pt z ydmdy Meats 4Dt RLBW thus Rt and the area of BRlttgt the form Rt x4aDt At 4th2 330 The boundary of BRm can be regarded as the front of population wave and that is where new neighborhoods and new shopping centers are being built The speed of the expanding of the metro area is R t x4aD 331 The analysis of the formation of the modern metropolitan has been carried out by Skellam Ske51 in early 1950 s for the spreading of muskrats Ondatra zibethz39ca in central Europe Today the muskrat is quite common in Europe and Asia lts range extends from Sweden and France in the west to the major river systems of Siberia in the east But muskrats are a recent addition from the New World According to a study published in 1930 by Ulbrich referenced by Skellam Ske51 in 1905 several muskrats found their way to freedom in a wilderness near the Moldau River about 50 kilometers southwest of Prague Radial dispersion and exponential growth followed Although the muskrat has many natural predators moat notably the mink in this instance there was evidently no natural imposed carrying capacity during the ensuing years the muskrat population grew rapidly and dispersed widely Over the years Ulbrich kept records of the spread of the muskrat His map of the locations of the dispersion front is shown 7 Cumulative areas de ned by the level curves contours of the map are listed in the following table 34 DIFFUSIVE MALTHUS EQUATION AND BIOLOGICAL DISPERSION 43 Year 19051909 1911 1915 1920 1927 Area ka 0 5400 14000 37700 79300 201600 We use a Maple program to t the above data into a function At ktz gtrestartwithstats withplots Times O 4 6 10 15 22 gtAreaz 0 5400 14000 37700 79300 201600 gtData00 45400 614000 1037700 1579300 22201600 gteqfit fit leastsquare x y ykxquot2 k Times Area 119777300 2 56 39 y 296433 3 119777300 1 Thus aD 32154229407 and the speed of the expansion is R t x4aD 296433 1134093989kmyear 47139 The equation 325 is the simplest model for biological invasion7 and it gives good estimate of the speed of invasion wave in the initial stage of the dispersion However it is also unrealistic in long term because of exponential growth In the next Chapter we will discuss more realistic invasion models 1 to DJ F Chapter 3 Exercises Consider the diffusion equation 74 Dum Assume that 74157 z 0zxf derive an ordinary differential equation satis ed by v and show that the solution utz is given by 318 Consider the convectivediffusion equation 2 8071380 V814 57 W7 a tgt0m Ru0zM60m 332 Use a change of variable 2 z 7 Vt to show that 74157 2 satis es ut DuM7 and derive the solution formula of 332 M as 7 Vt2gt u hm em if l 4th p lt 4D 333 In the fuel spill example in Section 337 if we assume now the canal current is V 001ms from the spill point toward Prof Shi s house7 then ftz satis es 332 Use the results in problem 2 to recalculate the questions in Example 31 In this case7 the center of the fuel patch is moving with a velocity V In the fuel spill example in Section 337 we assume that the current is from Prof Shi s house toward the spill point Then it is possible that the fuel never reaches Prof Shi s house if the current is fast enough Here reach means the house is inside the 120 of the center of the patch Determine for which V that will happen I4 4 F N 1 of 30 Once Upon a Time Constant diffusion J DVu Bu VJDAu 6t But is this truly realistic I4 4 F bl 2 of 30 Any initial condition even a point distribution instantly quotspreads outquot to cover an infinite domain Consider the onedimensional case 21 D 5921 I4 4 F N 3 of so I4 4 F bl 4 of so How Can We Improve This Use a nonconstant diffusion term J Du V ux t 5331 v J v D v at u u This makes intuitive sense in an insect population for example we would expect very densely populated areas to diffuse outwards more quickly than sparsely populated areas I4 4 F N 5 of 30 Of course now we need to figure out how to deal with nonconstant diffusion in our solution I4 4 F N 6 of so A General Approach Rewrite our equation as Eli v D v 0 at u u We can consider this to be an example of a general class of functions of the form Gx t u x t xx uxt utt 0 I4 4 F N 7 of so A General Approach The oneparameter family of stretching functions szax zebr H S H m S I4 4 F bl 8 of so A General Approach The oneparameter family of stretching functions a i5 x Ezebr zEcu a b and c are constants E is a real parameter on some open interval that contains 1 Define G to be invariant if there exists a smooth function f 6 such that GR f H Hi i Hm Hf H5 2 f5 G05 13 u x t xx uxt utt I4 4 F N 9 of so A General Approach Assume G is invariant This gives us GR f H f E Gx t u GE x Eb r EC u f5 Gx t u 65 x e t EC u f5 0 65 x e t EC u 0 Because G is homogenous I4 4 F bl 10 of 30 A General Approach Differentiate with respect to 62 3G 3G 61an thb cch 0 x 3t Bu Set 6 l which we can do because we restrict E to a domain that contains 1 Clever people look at this and see that the transformation we want to use is u tcb rz x z 2 FE I4 4 F M 11 of 30 3G 62 6G 62 3G 62 ax bt cu 0 92 3x 92 BI 92 Bu 92 1 92 I4 4 F N 12 of 30 A General Approach What have we accomplished with all our fancy math Gx t u x t xx uxt utt 0 ltgt gz r r39 rquot 0 I4 4 F N 13 of 30 A Specific Example Recall the problem we re actually working on Bu t V DuV u 6 Du D0 0 39quot I4 4 F bl 14 of 30 A Specific Example Letting m 1 gives us Bu 3 D0 Bu 3t bx uo 3x And because we re lazy we ll assume 30 2 1 so 0 Eu 3 Bu 3t 3x 3x The problem is now 9 9 9 15339 395 u 5quot ux 0 6x I4 4 F N 15 of 30 A Specific Example Other assumptions Since no organisms are being born or dying we require for all t gt 0 fmuOc t dx 1 and lim ux t 0 Xaioo I4 4 F bl 16 ofso GR f 7 if if 77 H7 7E f5 G05 13 u x t xx uxt utt i5 x Ezebr zgcu 3x a I9 W13 6071395 52072a 3 u 3 We have invariance if c b20 2agtc2a b I4 4 F N 17 of 30 A Specific Example u tcb rz x z 5 With our invariance condition u t2 Wb rz x z 5 I4 4 F bl 18 of 30 A Specific Example Let s be clever fmuOc t dx 1 2 aimb fC t Imrtabdx l termb fmrzdz 1 Timeindependence requires b 3 61 Which simplifies the transformation to u I m 12 2 xfll3 I4 4 F N 19 of 30 A Specific Example uztil mz gt3rr rzr390 z x t 13 This equation can be integrated to give 3 r r39 Z r constant Take the constant to be zero the solution is I4 4 F bl 20 of 30 A Specific Example Use our conditions to clean it up lim ux t 0 Xaioo means that Mug u w 0 x gt A And fmuOc t dx 1 means that 9 13 5 Now switch eveiything back to original coordinates I4 4 F N 21 of 30 A Specific Example 1 A2 W3 x2 x lt A t13 ux t 6 0 x gt A t13 About time Let s take a 100k I4 4 F N 22 of 30 A Specific Example Pretty Pictures I4 4 F N 23 of 30 A Specific Example The key feature of this solution is the sharp wave front at x f A t13 This wave is moving with speed 311 1 A fez3 dr 3 I4 4 F N 24 of 30 Comparing Constant and DensityDependent Diffusion at t 1 I4 4 b N 25 of 30 Comparing Constant and DensityDependent Diffusion at t 10 I4 4 F bl 26 of 30 What About the NotSimple Case You Ask Recall that the general form is Bu t VDuVu a Du D0 0 39quot and we assumed m l for all the work we just did Is there a general solution I4 4 F bl 27 of 30 What About the NotSimple Case You Ask Yes and here it is ux r 1 395 z21m lxl 5 r0 m 0 x gt r0 Mt Where I 12m Mr 2 to 1 3 0 z QFra to z 151 n12 MO r71 l1 2D0m2 D0 and no are positive constants Q is the initial density at the origin and r0 comes from requiring that the integral over the domain at all times be equal to Q I4 4 F N 28 of 30 I4 4 F N 29 of 30 I4 4 F bl 30 ofso Pay No Attention to the Man Behind the Curtain In30 ux t 1 WhichAbsx lt tl3 t23 x2 True 0 6t 1 x In31 uConstantDx t Exp V4xt 4t In39 Plotux 1 ux 2 ux 3 x 1 1 Plotstyle gt RGBColor1 0 0 RGBColor0 1 0 RGBColor0 0 1 In40 Plotux 1 uConstantDx 1 x 1 1 Plotstyle gt RGBColor1 0 0 RGBColor0 1 0 RGBColor0 0 1 n46 n47 n48 n49 n50 Plotux 10 uConstantDx 10 x 10 10 PlotRange gtAll Plotstyle gt RGBColor1 0 0 RGBColor0 1 0 RGBColor0 0 1 t m1 t0m At Gammai 32 r0m m WGama1m1 r0m2m t0m 2 m 1 uGeneralx t m W39hichAbsx s r0m At m 1 x 2 lm 1 True 0 At m r0m At m Periodic solutions A periodic solution is a solution 1tyt of da d d gxy such that 1t I T 513t and ytT yt for any 25 where T is a fixed number which is a period of the solution Periodicity in biology life proceeds in a rhythmic and periodic style the periodic patterns are not easily disrupted or changed by a noisy and random environment First example dy Polar coordinates 1 rcose y rsin 0 r O 9 1 Second example d d d fyx1 x2 y2 d Z xy1 x2 y2 Polar coordinates r r1 r2 0 1 Logistic equation for periodic solutions General polar coordinate equation da d a Mon may d Z wrx Amy Polar coordinates r M0 0 wr Examples Predator prey models lynx hare population shark fish population nerve conduction Hodgkin Huxley equations Fitzhugh Nagumo model chemical reactions Lokta Volterra predator prey eq uation da b t a13 my dy d t ca my Critiques of the equation there are infinitely many periodic so lutions and all of them are neutrally stable Ideally the system should have a limit cycle which all solutions tend to Limit cycle a periodic solution which is stable if the initial point is inside and near the cycle then the solution starting from there spirals out and tends to it if the point is outside and near the cycle then the solution starting from there spirals out and tends to it Predator prey model with Holling type II predation saturating interaction du am d2 duv d8 u U ub7 d8 cvub 1 u nullcline u 011 Z M a v nullcline 901 u be d c Equilibrium points 00 10 dbc bddcbc c ad c2 Dynamics when a d 1 u nullcline u 0 v 1 uu b be 1 c Equilibrium points 00 10 v nullcline v 0 u be b1 c bc 1 c 1 c2 39 1 1 b and 10 is a sink which attracts any initial values Case 1 cgt Two equilibrium points at 00 and 10 1 Case 2 0 lt c lt 1 H Three equilibrium points at 00 10 uovo bCC7 b1 CC 2bC 1 Stability of 110220 A brief theory of periodic orbits Theorem 1 Inside a periodic orbit there is at least one equilib rium point Corollary 2 Let 0340 be an equilibrium point If there is a stable or unstable orbit of this equilibrium going to infinity as t gt oo then there is no periodic orbit around this equilibrium point Corollary 3 Let 0340 be an equilibrium point If the direction of vector field on the nullclines near 0340 is not clockwise or counterclockwise then there is no periodic orbit around this equilibrium point Theorem 4 Poincar Bendixon If there is a donut region a region without equilibrium point that the vector field is pointing inside at any point on the boundary then there is a periodic orbit in the donut region Corollary 5 If in a region there is only one equilibrium point which is a source or spiral source and the vector field is pointing inside at any point on the boundary then there is a periodic orbit around this equilibrium point A typical case the solutions from outside are spiraling in and the equilibrium point is a spiral source Theorem 6 Poincare Andronov Hopf bifurcation If a bifurca tion occurs such that an equilibrium point changes from a spiral sink to spiral source and the solutions from outside are always spiraling in then there is a periodic orbit around the spiral source Jules Henri Poincare 1854 1912 One of greatest mathematicians in 20th century httpwww gapdcsst andacuk historyMathematiciansPoincareh1 8 du uv do uv 1 d5 M u ub d5 cvub 2 Stability Of u0v0 C1 c bc b C J 1 c 1 c bc O c1 c bc bgt Eigenvalue equation A2 1 c 1 c bc 2 O 1 b S iraI sink or sink 1 b bltO or gt D c c C 1b 1 b Spiral source or source 1 c 96 b gt O or O lt c lt 1 b 9 Conclusion c is where Hopf bifurcation occurs When 1 gt 1b but near 1b is a spiral sink and all solu c u v 1b 1b 0 0 tions tend to 110220 in an oscillating fashion when c passes 1 b but near 1 b 1 b solutions away from 110220 still tends toward 110220 spiral in ward thus a periodic solution with small oscillation emerges around uovo and it is a limit cycle 110220 becomes a spiral source but all Examples b 05 the Hopf bifurcation point is c 13 m 033333 b 02 the Hopf bifurcation point is c 23 m 066667 10 Theorem proved around 1980 For ul u z cv I m d 7 7 s u I b d5 u I b 1 1 b if gt then 10 is as m toticall stable if lt lt 0 1 y I0 y 1 b 6 1 1 1 b then is as m toticall stable and if 0 lt c lt 11 U0 110 y D V 11 then there exists a unique periodic solution and it attracts all positive initial values except the equilibrium points 00 10 and u0 00 Hilbert problem 16 Consider u fuv 12 guv If f and g are polynomials prove the system has only finite many limit cycles When f and g are quadratic functions prove there are at most three limit cycles Counterexample is found by S Shi in 1979 four limit cycles 11 Turing instability i iluum u1 71 Jib 1 6 01 m vvm cvu b 5136 01 4 ux07t u17t vx07t 11361 2 0 uaO u0a Z O vaO 12033 Z O 13 6 01 Nope since 91 O we need fugv lt O for Turing instability 1 b But many Turing type bifurcations occurs when 0 lt c lt 1 b along the equilibrium solution u0v0 The global dynamics of the reaction diffusion system is still not known 12 Neural modeling excitable systems FitzHugh Nagumo system evv al v w dt dw b 5 v c w dt 2225 excitation variable wt potassium conductance a 01 b 05 c 0 e 001 equilibrium 00 is asymptoti cally stable but it is an excitable system where a large excursion occurs if a perturbation is beyond a threshold a 01 b 05 c 01 e 001 equilibrium 010 is unstable and a limit cycle exists relaxation oscillation 13 v TIM Iollt gr u VVILLIAM MARY Pattern qumation in Chemical and Blologleal Systems JUNPING SHI E iP College of William and Mary Williamsburg Virginia 231 87 CSUMS Lecture College of William and Mary February 18th 2009 Bi ircation in 11 ReactionDiffusion Equations Many biological processes can be simulated by mathematical models involving temporal and spatial variables Example a video of cell division a video of simulation of a reactiondiffusion system Reactiondi usion systems are mathematical models that describe how the concentration of one or more substances distributed in space changes under the in uence of two processes local chemical reactions in which the substances are converted into each other and diffusion which causes the substances to spread out in space Bifurcation 7 D 21 Mathematical Models W DlAuJct fuxtv7ta W D2A21Jct 9uatavxvt 15 time variable x x1 x2 spatial variable ux t UZE t density of substances at time t and location 16 82u1ct 82UZE t 996 9963 Auxt divVuEt Au Av are Diffusion transport of molecules from a region of higher concentration to one of lower concentration by random molecular motion f u v gu v are Reaction deathbirth chemical reactiongeneration Bifurcation 7 D 31 Alan Turing 19121954 39 One of greatest scientists in 20th century 39 Designer of Turing machine a theoretical computer in 193 Os 39 Designing electromechanical machine Which breaks German Uboat Enigma helping the battle of the Atlantic 39 Initiate nonlinear theory of biological growth Turing 1952 The Chemical Basis of Morphogenesis Philosophical transaction Royal Society of London Series B 237 httpwwwturingorguk Bifurcation n 4 l Turing s idea ODE 12 u aw 11 9a 11 Reactiondiffusion system 2 at d1 Au fu 21 ct dgAU gu 21 Here ux t and UZE t are the density functions of two chemicals morphogen or species which interact or react 39 A constant solution ut x uo 2125 x 210 can be a stable solution of l but an unstable solution of 2 Thus the instability is induced by diffusion 39 On the other hand there must be stable nonconstant equilibrium solutions or stable nonequilibrium behavior which have more complicated spatialtemporal structure httpenwikipediaorgwikiMorphogen Bifurcation 7 D 51 Oscillatory Chemical Reactions An example of oscillatory chemical reaction don t do this in home though Chemical reactions were believed to always reach equilibrium states even when the reactions are reversible But the discoveries since 1960s con rm the eXistence of chemical oscillations It is now believed that every living system contains hundreds of chemical oscillators Bifurcahun 7 n 61 Equilibrium solutions and Periodic solutions For a dynamical system ODE reactiondiffusion model etc are the ones independent of time Periodic solutions are the ones which repeats itself in time They usually dominate the long time behavior of the dynamical system l l I W i l l l l quot My M aili KL L n n n in v E a Bifurcation 7 p 71 Hopf Bifurcation Theorem Consider ODE x fA7 at A E R x E R and f is smooth i Suppose that for A near A0 the system has a family of equiliina x0 if Assume that its Jacobian matIix AA fw A7 x A has one pair of complex eigenvalues MA i iwA MAO 0 MAO gt 0 and all other eigenvalues of AA have nonzero real paIts for all A near A0 Ifu A0 74 0 then the system has a family ofpeIiodic solutions As7 for s E 07 6 With peIiod Ts such that As a A0 Ts a 27rwAo and 7 x0AoH a 0 as s a 0 laxmm 7 n 21 LotkaVolterra model l 1 Alfred Lotka 1880 1949 Vito Volterra 1860 1940 Cell q ua bu cm CL dv fuv Q35 Bifurcation 7 D 91 Functional response d ZL ua bu c uv v a dv f uv u predator functional response u u LotkaVolterra u 1 umu Holling type H m the handling time of prey Holling 1959MichaelisMent0n biochemical kinetics Biological work RosenzweigMacArthur American Naturalist 1963 Rosenzweig Science 1971 Paradox of enrichment May Science 1972 Existence and uniqueness of limit cycle Bifurcation 7 D 101 Basic analysis of the model d u u1 u muv dv muv dt au dt au 1 Nullclineisocline u 0 v w v 0 d mu m au d Solvingd u onehaveuAE a au m d 1 A A Equilibrium pomts 00 10 Ana where m W We take A as a bifurcation parameter Case 1 A 2 1 10 is globally asymptotically stable Case 2 1 a2 lt A lt1 10 is a saddle and Ana is a locally stable equilibrium Case 3 0 lt A lt 1 a2 10 is a saddle and A m is a locally unstable equilibrium A 1 a 2 is a Hopf bifurcation point Bifurcation 7 D lll m t 1 W lt A lt 1 1m g 9 km g all tax lt A lt 1 W2 LED g am39aro colecag gmd 2A9 mix a ma y um b b qwi iibzr mm mama a Emit 11m 7 w H311 Summary duu1umuv vmuv dt au dt au 1 Nullclineisocline u 0 v w v 0 d mu m au d Solvingd u onehaveuAE a au m d 1 A A Equilibrium pomts 00 10 Ana where m W We take A as a bifurcation parameter M A 2 1 10 is globally asymptotically stable M 1 a2 lt A lt 1 A m is a globally asymptotically stable 0 lt A lt 1 a 2 the unique limit cycle is globally asymptotically stable A 1 a 2 is a Hopf bifurcation point Bifurcation 7 D 1 31 New result of this ODE Hsu Shi 2009 Hsu SzeBi Shi Junping Relaxation oscillator pro le of limit cycle in predatorprey system Disc Cont Dyna SystB Motivated by numerical observation Graph of limit cycle Parameters a 05 m 1 d 01 118 x 0056 period T m 37 A review of stability and dynamical behaviors of differential equations scalar ODE at fu system of ODEs at u v vi 9uav7 reaction diffusion equation at DAu fu 1 E Q with boundary condition reaction diffusion system D A ut u u v 1 E Q with boundary condition vi DyA U gu 1 All equation is in form of Ut FU where u can be a scalar or vector spatial independent or dependent Abstract Equation Ut FU Equilibrium solution U0 such that FUO O linearized operator F UO U0 is stable if the eigenvalues of equation F U0w Aw are all with negative real parts Linear behavior Since the equation Ut FU is approximately the linearized equation Ut F UU U0 near the equilibrium solution U 2 U0 then U05 m ZCiexpOwtkbi near U 2 U0 where Moi are the eigenvalue eigenvector pairs of F U0w Aw If A1 max M then U05 m 01 expgt1tglgt1 if UO m U0 Turing bifurcations in reaction diffusion system at um Afuv vi dvm Aguv ital77quot 0 vt7 vt77r 0 If u0v0 is an equilibrium so that fu0vo gu0v0 0 then 110220 is also an equilibrium of at fuv vt guv Linearized operator ODE Jacobian J lt f f 9a 911 I A O fu f1 PDE diagduAadUA J lt o my gt lt 9a 91 Ash Afuuo vow Af ne vow Mb Eigenvalue problem dAgb Aguu0 vow Ag uo vow 2 Mb Chasm 3007 2 0 M140 2 15x07 2 O 3 If all eigenvalues of J are with negative real parts so u0v0 is a stable equilibrium solution for ODE is u0vo also a stable equilibrium solution for PDE It seems so since the additional part is consist of diffusion oper ators only and diffusion is supposed to stabilizing But as Alan Turing pointed out 110220 could be an unstable equilibrium solution for PDE even if it is stable for ODE So diffusion has an unstable effect for such system How is that possible Let s calculate now Ideas Eigenfunction could be in form of Vcosmx where COSk13 is the eigenfunction of a ugh O gb7r O and V is a vector Substituting my 2 VCOSk13 into the equation we get AJ k2DV uV fu f1 10 whereJ andDlt gu 91 0 d of the reaction diffusion system are ulyk uzk which are eigen values of J MD with k 12 Thus the eigenvalues Question Can some of ulyk uzk be positive Which one is the largest Classification of linear system Y lt gtY U U Two real eigenvalues 1 A1 gt A2 gt 0 source 2 A1 gt A2 O degenerate source 3 A1 gt O gt A2 saddle 4 A1 O gt A2 degenerate sink 5 O gt A1 gt A2 sink Two complex eigenvalues ME 2 ad bi 1 a gt O spiral source 2 a 0 center 3 a lt O spiral sink One real eigenvalue A1 A2 A 1 A gt 0 star source or trying to spiral source 2 A 0 parallel lines 3 A lt 0 star sink or trying to spiral sink Generic Cases most likely not fragile and nonlinear system behaves similar to linear system Source Sink Saddle Spiral source Spiral sink Stable sink or spiral sink real part of eigenvalues is negative Unstable source spiral source saddle at least one of eigen value has positive real part Trace determinant criterion for stability AQ AABO A2 4B gt 0 two real eigenvalues 1 A gt O B gt 0 source 2 A lt O B gt O sink 3 B lt O saddle A2 4B lt 0 two complex eigenvalues 1 A gt O spiral source 2 A lt O spiral sink Stable A lt O and B gt O Unstable A lt O and B gt O or B lt O Calculations 2 p Afu kQ gtfv AMI J k D lt gu MAgv k2dgt Equation of eigenvalues 2 Afugv k21du ange rm k2dfugvxk4d O or M2 Au B 0 where A Afu A91 k2 k2d B A2ltfugv fugu k2dfu 910A k4d We want it to be unstable thus either A lt O and B gt O or B lt 0 but we also want it is stable wrt ODE thus A fu 91 lt O and B fuqu fvgu gt O and it implies A gt 0 So we must have B lt O saddle type and only one of pm or u27k is positive 9 So we hope B A2fugv fvgu k2dfugv k4d lt o for some integer k Since fugq fvgu gt 0 then dfu 91 gt 0 otherwise B gt O for all Ak So we get necessary conditions of Turing instability fugv fvgu gt 0 fu l39gv lt 0 dfu l39gv gt O fu and 91 must be of different sign here we assume that fu lt O and gvgtO so Oltdlt 1 Now from A2fugv fugu 1132de 91 k4d lt 0 we solve d Agvk2 mew d lt 022 Ara where detJ quv fugu So Turing instability will occur if fugv fvgugt0a fUgUltOa degUgtOa fult07 911gtOa Oltgtxlt 2 2 gvk 7 and O lt d lt Agvk ADetJ DetJ k2k2 Afu When these conditions are met then k th mode is unstable COSk13 is an eigenfunction The band of unstable modes is given by k2 e a1a2 where a1 and a2 are the roots of equation A2ltfugv fugu 0de 910A da2 01 a 21d ltde 91 i de 902 4dDetJ 10 An artificial example J fu f1 2 3 2 gu 91 4 2 For ODE it is stable since A fu I gq 1 lt O and B Wealso havefu2 3ltOandgq2gt0 Hence the k th mode is unstable if 2 OltAltk2 Oltdltdkgt2 2A1 A 2A4 A 2A9 A d1A W dQW 44 3A d3m 99 3A 11 Parameter Ad regions for Turing instability Horizontal axis is A vertical one is d and the graphs are d1gt d2gt and d3gt 12 For example if Ad 201 then the band of unstable mode is 25431457 so k mode k 2345 are unstable From calculation of Maple we find that the reaction diffusion system has 4 positive eigenvalues The largest eigenvalue is M173 1114 with eigenvector V 03161273 Hence near the equilibrium point if we perturb the system we have 1605 N up 111415 0316 220 gt N lt 220 6 C0563 1273 A spatial pattern with period 27r3 is generated by Turing in stability The characteristic wave number is k 3 and the characteristic wave length is 27r3 13 Turing bifurcations When the parameters Ad change then k mode can turn from stable to unstable For example in our example when A 2 and d decreases from 03 to 0 it cuts through d200 d3gt d4gt and d5gt Each time it cuts a dkgt k mode becomes unstable Similar phenomenon occurs if we changes A Each time a k mode becomes unstable a curve of equilibrium solutions emerges from the constant equilibrium solution The properties of these curves are still unclear for most reaction diffusion systems People assume that the system will generate a spatial pattern with the most unstable mode of the constant solution but this is not verified mathematically However the Turing instability helps us when we do numerical simulation of the system 14 Chapter 1 Derivation of ReactionDiffusion Equations 11 Fick s Law Diffusion mechanism models the movement of many individuals in an environment or media The individuals can be very small such as basic particles in physics bacteria molecules or cells or very large objects such as animals plants or certain kind of events like epidemics or rumors The particles reside in a region which we call 9 and we assume that Q is open subset of R the n th dimensional space with Cartesian coordinate system with n 2 1 In particular we are interested in the cases of n 1 2 and 3 but most material here are true regardless of the dimensions of the space sometimes 71 1 and n 2 2 may be different as we will see The main mathematical variable we consider here is the density function ofthe particles Pt m where t is the time and z E Q is the location The dimension of the population density usually is number of particles or organisms per unit area if n 2 or unit volume if n 3 For example the human population density is often expressed in number of people per square kilometer A list of world population and population density can be found at http en2 wikipedia orgwikiListofcountriesbypopulationdensity However in such data table we can only nd the population density for countries and a country is not a point on our map In reality population density is always associated with a scale like country city county town and street But as in many other mathematical models we will assume that the function Pt x has nicer properties like continuity and differentiability which is in fact reasonable when a population with a large number of organisms is considered Technically we de ne the population density function Pt x as follows let x be a point in the habitat Q and let On 301 be a sequence of spatial regions which have the same dimension as Q surrounding m here On is chosen in a way that the spatial measurement Onl of On length area volume or 1 2 CHAPTER 1 DERIVATION mathematically the Lesbgue measure tends to zero as n a 00 and On 3 On then Number of organisms in On at time t 11 KM if the limit exists Realistically as long as the scale of data collecting is small enough the population density function is usually very well de ned Clearly the total population in any sub region O of Q at time t is Pltmgt 322 Ptxdz 1392 0 The question we are interested now is how the function Pt z changes as time t evolves and as the location z varies Population can change in two ways one is that the individual particles can move around and the second is they may produce new individuals or kill existing individuals due to physical chemical or biological reasons We shall model these two different phenomena separately How do the particles move In general that is a highly complicated process which can be attributed to a lot of reasons For example the reasons of the emigration of human can be looking for a better life looking for a better job political reason or religious reason etc Let s review a few historical emigration trends since the sixteenth century European emigrated westward to America in the last two centuries Asian emigrated eastward to America inside United States people move to west from east and people move to Rocky mountain area from eastern or western coasts Although economical reason can be quoted as motivation to move we can see clearly in all these cases people move from areas where population density is high to areas where it is lower This is similar to many physical phenomena like the heat transfer from warmer place to colder place or the dilution of chemical in water A Chinese proverb is People goes to high place water ows to low place It is a natural phenomenon that a substance goes from high density regions to low density regions The movement of Pt z is called the ux of the population density which is a vector The high to low77 principle now means that the ux always points to the most rapid decreasing direction of Pt m which is the negative gradient of Pt This principle is called Fick7s law and it can be represented as Jt z 7dxVmPtz 13 where J is the ux of P dz is called diffusion coe icient at m and Wm is the gradient operator On the other hand the number of particles at any point may change because of other reasons like birth death hunting or chemical reactions We assume that the rate of change of the density function due to these reasons is ft x P which we usually call the reaction rate Now we derive a differential equation using the balanced law We choose any region 0 then the total population in O is f0 Pt zdm and the rate of change of the total population is d P t d 14 d O lt 790 z lt gt The net growth of the population inside the region 0 is m zPtmdz 15 O 12 RANDOM WALK 3 and the total out ux is Jt x nzdS 16 30 where 80 is the boundary of O and nz is the outer normal direction at m Then the balance law implies d Ptmdm 7 Jtm nxdS ftmPtxdz 17 dt 0 30 0 From the Divergence Theorem in multi Variable calculus we have JtxnmdS divt 18 30 O Combining 13 17 and 18 and interchanging the order of differentiation and integration we obtain Wm dwdzVmPtz ftzPtx dz 19 0 at 0 Since the choice of the region 0 is arbitrary then the differential equation W divdxV1Ptx ftxPtx 110 holds for any 25 The equation 110 is called a reaction diffusion equation Here divdzVmPt is the diffusion term which describes the movement of the indiViduals and ftzPt is the reaction term which describes the birth death or reaction occurring inside the habitat or reactor The diffusion coef cient dz is not a constant in general since the enVironment is usually heterogeneous But when the region of the diffusion is approximately homogeneous we can assume that dz E 1 then 110 can be simpli ed to 88 dAPftzP 111 11 82P where AP divVP E W is the Laplacian operator When there is no reaction occurs the 1 1 equation is diffusion equation 8P dAP 112 a lt gt In classical mathematical physics the equation Tt AT is called heat equation where T is the temperature function So sometimes 111 is also called a nonlinear heat equation Conduction of heat can be considered as a form of diffusion of heat 12 Random Walk A random walk considers a walker a particle or a rabbit which starts at a point and takes steps in a random direction Sometimes the steps can also be of random length as well The random walk can take place in a plane along a line or in higher dimensions 4 CHAPTER 1 DERIVATION The simplest random walk considers a walker that takes steps of length Am to the left or right along a line and after each At time units the walker will take one step If the walker is at location m0 at the time to then at time t to 1 At the walker will either be at me 7 Am or me 1 Am Normally the chances or going left or right should be equal thus the probability of the walker going left or right is 12 Now we assume that many walkers are walking with the same time frame simultaneously with the same step size and on the same lattice on the line We de ne Pt m as the number of walkers at time t and location m Then after one time step At everyone who is at me is gone now go to me 7 Am or me 1 Am and half of those who are at m0 7 A and half of those who are at me A move to m me now So we have 1 1 Pt0 1 At 0 Pt00 7 Pt0 0 We use the Taylor expansion of the function Pt m 8P 182P Pt0 1 At 0 Pt00 Et00At Wt00At2 1 1 1 18p 182P 2 Pt00 7 Am 7 Pt00 Ea o 07A ZWI 00A 115 1 1 18p 182P 2 Pt00 Pt00 55050 imam 1 By substituting 115 116 and 114 into 113 we obtain 8P 1 82F t At t A 2 117 at 07950 2 atg 07950 95 t Here we assume that both At and Am are small quantities thus the higher order terms in the Taylor expansions the parts are smaller terms compared to the two remaining terms in 117 By dividing At we now have 8P Am2 82F 2 at t00 2At atz t00A 118 Now we assume that A902 2At The last assumption can be satis ed by designing the scale of the random walk and when the limits in 119 are taken we arrive at A2570 A3570 and HDgt0 119 8 8 8 t070 DW7 0707 13920 P 2P t the one dimensional diffusion equation From the above argument the diffusion constant d have the dimension as Am2At or distance2 121 time 13 REACTION 5 The diffusion constant can often be estimated through experiment since 121 implies the time that the substance diffuses is proportional to the square of the distance that the substance diffuses 239 e T N D2 122 where T is the time that the substance diffuses and D is the distance that the substance diffuses We will get a more precise conclusion later with calculations but we can see from the following chart of the diffusion constant of oxygen that D depends on the nature of the substance the media where it is diffusing and the temperature Temperature 00 Media D cmZsec 0 air 178 X 10 1 20 air 201 x 10 1 18 water 241 x 10 5 25 water 458 x 10 5 The assumption that the random walk is unbiased may not be true for certain cases and when one direction in the random walk is favored an equation with a rst order derivative is the result see homework problem 8P 8P 82F V D 123 875 1 8m 8752 l The rst order derivative term is called a convection term and 123 is a diffusion equation with convection In higher dimensional space the random walk can be biased on any direction and the general reaction diffusion equation with convection is 814 8214 EVVu7DWftxu 124 where 397 E R is a vector indicating the direction of the biased motion see homework problem Such equation is more realistic when the media of the diffusion moves to a directionifor example the wind or the water ow 13 Reaction In the reaction diffusion equation 8P at dAP m x P 125 the ft m P represents the birthdeath or reaction process as will see in the following chapters reaction and diffusion both contribute to the interesting dynamical behavior of the solutions of the equation and sometimes it is intriguing to compare the dynamics of 125 with its more familiar and less sophiscated cousin 8P 5 7 m P 126 6 CHAPTER 1 DERIVATION where P Pt is usually not the density function but the function of total number of the particles at time t Generic population model usually assume that fP kP Malthus linear growth fP kP lt17 Logistic growth 127 The analytic methods for 126 with these growth rate functions are well known in the undergrad uate differential equation textbooks Reaction diffusion equations corresponding to these reaction rates are 8P E DAP kP Diffusive Malthus equation 128 3P P E DAP kP lt17 Diffusive logistic equation 129 0 4 l 0 3 1 0 8 0 2 0 6 01 04 W 0 2 0 4 u0 6 0 8 0 2 r 1 70 2 T 0 2 0 4 u0 6 0 8 1 02 04 06 08 Figure 11 a Logistic upper b Weak Allee effect middle growth rate left growth rate per capita right Similar to the non spatial versions 128 and 129 are equations describing the evolution of a spatially distributed population which satis es generical growth patterns In these two examples if we rewrite the equation as 2 1 DAP P903 130 then the growth rate per capita gP are constant 91P k for Malthus equation and decreasing linear function 91 P k1 7 PN for logistic equation 13 REACTION 7 Figure 12 Srong Allee effect growth rate left growth rate per capita right In population biology for some species a small or sparse population may not be favorable since for example mating may be dif cult This is called Allee effect Mathematically the growth rate per capita function gP will not have the maximum value at P 0 If gP is negative when P is small we call such a growth pattern has a strong Allee effect A typical example is 412 17 lt71gtP9P 131 where 0 lt M lt N M is the sparsity constant and N is the carrying capacity If the growth rate per capita gP is smaller than the maximum but still positive for small P we call such a growth pattern has a weak Allee effect A typical example is dP P P W ikP lt17Ngt M71 7PgP 132 where M lt 0 lt N The corresponding reaction diffusion equation with Allee effect is 8P P P E DAPaPlt17Ngt M71 133 whereagt0Ngt0andNgtM Chemical mixing and reaction can also be modeled by reaction diffusion equations We rst consider a typical mixing problem77 in ordinary differential equation course A tank currently holds 100 gallons of pure water A solution containing 2 lb of salt per gallon enters the tank at a rate of3 gallonminute Assume that the solution is well stirred so the salt is evenly resolved in the solution A drain is opened at the bottom of the tank so that the volume of the solution in the tank remains constant How much salt is in the tank after 60 minutes How about after a long long time Let At be the concenration of the salt the tank Then by the balance law At satis es an ordinary differential equation dA E 006 7 0031 A0 0 134 8 CHAPTER 1 DERIVATION The solution of the equation is At 217 e gtloo and as t a 00 At a 2lbgallon While the model 134 correctly describes the mixing as an exponential approach to an equilibrium state it is a simpli cation of the real chemical processiespecially the spatial distribution of the chemical is ignored and the chemical is assumed to be evenly resolved in the solution immediately which is also not true We use the reaction diffusion equation to model the problem again The tank which holds 100 gallon is the spatial domain 9 in R3 which can be a cylinder or cube shaped the concentration of the salt in the water is At m where t 2 O z E 9 we assume that the salt molecules moves randomly in the media water Since there is no reaction involved here the evolution of the density function A can be described by a diffusion equation with an initial condition DAA A0z 0 135 Note that 135 does not include the addition from the in owing solution which contains 2 lb per gallon of salt and the subtraction from the drain In fact these additions or subtractions occur not in the domain but on the boundary of the domain For example the in owing solution comes into the tank from the solution surface and the drain is on the bottom of the tank We shall discuss such boundary conditions of the problem in the next section In the above problem if salt is replaced by a chemical which will react say the name of chemical is A and the chemical reaction is Adv 01 136 where C and D are two other kinds of molecules and k is the reaction rate the percentage of B molecules which will react in a unit time and a unit volume Then the reaction diffusion equation is 8A DAA 7 kA 137 8t where At z is the concentration function of A Of course most chemical reactions involve more than one chemicals thus a system of reaction diffusion equation is required to describe such reac tions For example we consider a chemical reaction A B J 2A 138 This is called an autocatalytic reaction as one A molecule acting as catalyst Let Atm and Bt x be the concentration of A and B respectively Since over all one B molecule is consumed and one A molecule is generated in the reaction the equations are 8A 8B DAAAkAB DBABikAB 139 8t 8t Here we assume the possibility of A and B molecules colliding in a unit volume of media is AB and reaction rate k is de ned as the percentage of reaction when a collision occurs The possibility of a multiple molecule reaction occurring is proportional to the multiplication of concentration of each reactant For example consider pA B J p 1A 140 14 BOUNDARY 9 where p gt 1 is an integer Then similar to 139 we have DAAA MW 88 DBAB 7 kAPB 141 A simple reaction which has been studied extensively is Gray Scott model Consider reactions A 213 711133 B 7 12 C 142 and we assume that A is continuously and evenly feeded into the reactor and A B and C are continuously and evenly removed from the reactor all with the same rate k3 Then the reaction diffusion model is 8A 82A EDAW 2 7 4132 kg 7 kgA 88 D3887 klABZ 7 sz 7 kgB 143 We will come back to 143 in Chapter 4 to discuss the behavior of the solutions 14 Boundary Conditions As we can see from the example of chemical mixing model the boundary conditions are as impor tant as the differential equations in the model In the situation of a chemical reaction additions and subtractions of chemicals can be included in the equation but often they are imbedded in the boundary conditions The simplest boundary condition is when there is no additions and subtrac tions of any chemicals through the boundary and all chemicals are generated inside the reactor and they remain there for example the chemical reactions occurring in a Petri dish with an im permeable wall Such system is called a closed one For each chemical involved if Jt z is the ux of the chemical then the ux across on a boundary point z is Jt z nm and if we assume the Fick s law 13 then for a close system at each boundary point m Vut m nz 0 144 where utz is the concentration function of the chemical Therefore a well posed which means solution is uniquely determined under reasonable assumptions closed reaction diffusion equation is an initial value boundary value problem 814 a dAu ftmu tgt 0 m E 9 Mom 7 W m e n 145 Vutm nz 0 tgt 0 m E 89 Here 89 is the common notation used for the boundary of 9 Closed reaction diffusion systems are similar In mathematical physics or uid dynamics 144 is called no ux or homogeneous Neumann boundary condition The second one is named after Germany mathematician and physicist Franz Ernst Neumann 1798 1895 The closed system is also relevant in many population ecological models For example an isolated island can be considered as the spatial domain of the system if all living species on the island do not attempt to emigrate from the island then the 10 CHAPTER 1 DERIVATION ecosystem of the island can be considered as a closed one For the classical heat conduction equation Tt dAT the no ux boundary condition means the boundary is heat insulated The no ux boundary condition can also be called as re ecting since particles diffusing to the boundary can be considered as re ected back In the case of spatial dimension n 1 the domain 9 is an interval a b then the no ux boundary condition becomes 814 814 05612 021 0 146 Now we can come back to the chemical mixing model considered in Section 13 The system there is not closed We assume that the in ow comes in from a portion F1 of the boundary 89 and the out ow comes out from another portion P2 Then the total ux across F1 in a unit time is J2522 nxds 6 147 F1 and we assume that the in ow is evenly distributed on P1 then for each x 6 P1 6 Jtm nz DVut m nz W 148 1 where D gt 0 is the diffusion constant and P1 is the surface area of P1 Thus 6 VA2522 2222 m 25gt 0 22 6 r2 149 Similarly 3A 25 VAtz 2222 7 13 25gt o 22 6 r2 150 For the other parts of the boundary the no ux boundary condition is appropriate Therefore the full initial value boundary value problem is 8A 82A majm tgt0xe 2 u0z 0 m E Q VA2522 1222 25gt 0 22 6 n 151 D P1 7 VA25 25 nz A25z 0 25gt 0 z 6 P2 Dngl VAtx nz 0 tgt 0 m E XIPl U P2 The boundary condition on F2 is a special case of homogeneous Robin boundary condition In general there are three commonly used boundary conditions NeumannVutz nz gtm tgt 0 z E 89 152 Robin Vut m nz axut m tgt 0 m E 89 153 and Dirichlet utz gtm tgt O z E 89 154 14 BOUNDARY 11 In the Robin boundary condition az 2 O and gtz is a function de ned on the boundary 89 As we can see from the previous example gtz may not be a continuous function so we assume that j is an integrable function The precise and optimal setting for such problems requires much more mathematical knowledge thus we will not go deeper in that direction Often we consider the case when gtz E O the homogeneous boundary condition When gtz is not zero then the boundary condition is nonhomogeneous Again as we can learn from the previous example it is unrealistic that the boundary condition is uniformly one type on the whole domain 89 In fact the boundary conditions in 151 is a mixed boundary condition which is one of the three boundary conditions on each of several components of 80 Dirichlet boundary condition assumes that the solution presumes the value gtz for each point on the boundary For many cases the homogeneous Dirichlet boundary condition makes the math ematical analysis easier so in many examples of this book we will consider that case For most chemical reaction no ux is more appropriate than Dirichlet boundary condition When the prob lem assumes Dirichlet boundary condition ut z E c a constant we can think ut z takes the value c for all the points outside of 9 For example for the heat equation Tt DAT constant Dirichlet condition means the outside environment has a constant temperature To For ecological applications sometimes we can assume homogeneous Dirichlet boundary condition ut z E 0 on 89 which implies that any individual of the species cannot survive outside of Q on even on the boundary of 9 Thus in this case the boundary is lethal sometimes called absorbing and any individual who wanders outside due to diffusion is killed by the sterile exterior environment or the deadly boundary A good example perhaps is the high voltage electric fence around the Jurassic Park When studying the periodic patterns generated by the reaction diffusion mechanism we will encounter periodic boundary condition First of all for such problem the domain and the media of diffusion both should have a periodic structure For example the onedimensional diffusion equation with periodic boundary condition is 814 8214 D t gt 0 E b 8 832 m a gt 14071 We 90 6 M uta utb tgt 0 m 6 ab u u 0201 tb tgt 0 m 6 ab 155 When the function u satis es 155 we can extend u spatially to z E 700 00 periodically The extended function satis es the differential equation for all z E 700 00 since the function values and the derivative values at the boundary match and the second derivative values will also match because of the equation Finally in some models of physics and biology we may consider the reaction diffusion equation on the whole space R The boundary of the whole space is z 00 or more precisely the boundary condition is the limiting behavior of the solution when a 00 where is the the distance from x to the origin In that situation a natural condition is that lim utm 0 l 1 Waco l l im Vutm 0 156 12 CHAPTER 1 DERIVATION It means that the chance of the dispersal of the articles to in nity is zero H to 03 F 01 Chapter 1 Exercises Let a fm y sinmcosy and b fm y 3x2 5342 a Calculate Vf b Calculate A c Use Maple to draw the graphs of fm y and the level curves of fz y contour graphs Let fm yz coszyezz a Calculate Vf b Calculate A c Show that Af2 2V Vf 2fAf and use that to calculate Af2 From satellite pictures the forest service nd that a new species of wild grass is spreading fast A data matching turns out that the density function Pt m y is approximately 1 2 yz Pt my Re t 475 157 a Verif that P t z satis es a reaction di usion e uation y 7 7y q 8P 82F 82F E W 1 8 242 2P 158 b Use the simulation function in Maple to simulate the spread of the this wild grass The two dimensional diffusion equation is 2 2 2 1 D 159 where D gt O Show that via a change of variables x rcos0 y rsin0 160 159 becomes 2 2 1 161 Let Pt x be the population density of species A at the location z and time t We assume that species A lives on a line and the movement of the individuals follows a random walk Suppose that in the random walk there is a bias to the right that is the probability of moving right is a gt 05 and the probability of moving left is b lt 05 but a 7 b is a small number 14 BOUNDARY a 004 to 13 Show that Ptz satis es a diffusion equation here we assume that there is no reaction with a convection term 162 where V and D are constants which you should de ne a Let Ptzy be the population density of species A at the location my and time t We assume that species A lives on a two dimensional space7 and the movement of the individuals follows a random walk Suppose that in the random walk7 there is no bias to any direction Set up a rectangular grid on R2 to model the random walk7 and show that Pt7 m y satis es a diffusion equation 8PD 811813 ati 8x2 8342 where D is a constant which you should de ne 163 Suppose that X7Y is another coordinate system in R2 which is a rotation of my X am 1 by7 Y ibm ay where a2 62 1 Show that Axyu Amyu This shows that the derivation of 163 is independent of the choice of the the coordinate system A T V c A more general diffusion equation is the anisotropic diffusion equation 8P 82F 82F E Dlw Dza yz 164 How do you change the assumptions in a to derive 164 Note that 163 is called isotropic diffusion equation in comparison with 164 Use the approach in Section 11 to derive the 124 The rst chemical respectable model is proposed by Prigogine and Lefever in 19687 which is later called Brusselator since Prigogine s lab was located in Brussels at that time The reactions are A 4 X B Y gtk2 Y D 2X Y 4 3X X 4 E 165 Suppose that A and B are feeded so that they keep constant concentrations in the reactions Write a system of reaction di usion equations for X and Y A tube with length L connects two reservoirs which contain salt waters The cross section yz direction of the tube is so small so we can assume that the concentration of the salt water is same for any point on a cross section Let the concentration of the salt water in the tube be Ct7 m 0 lt z lt L We assume that the diffusion of the salt is onedimensional in the tube 2 4 D2 3t amz Determine the boundary and initial conditions according to the following setups 166 Periodic solutions A periodic solution is a solution 1tyt of da d d gxy such that 1t I T 513t and ytT yt for any 25 where T is a fixed number which is a period of the solution Periodicity in biology life proceeds in a rhythmic and periodic style the periodic patterns are not easily disrupted or changed by a noisy and random environment First example dy Polar coordinates 1 rcose y rsin 0 r O 9 1 Second example d d d fyx1 x2 y2 d Z xy1 x2 y2 Polar coordinates r r1 r2 0 1 Logistic equation for periodic solutions General polar coordinate equation da d a Mon may d Z wrx Amy Polar coordinates r M0 0 wr Examples Predator prey models lynx hare population shark fish population nerve conduction Hodgkin Huxley equations Fitzhugh Nagumo model chemical reactions Lokta Volterra predator prey eq uation da b t a13 my dy d t ca my Critiques of the equation there are infinitely many periodic so lutions and all of them are neutrally stable Ideally the system should have a limit cycle which all solutions tend to Limit cycle a periodic solution which is stable if the initial point is inside and near the cycle then the solution starting from there spirals out and tends to it if the point is outside and near the cycle then the solution starting from there spirals out and tends to it Predator prey model with Holling type II predation saturating interaction du am d2 duv d8 u U ub7 d8 cvub 1 u nullcline u 011 Z M a v nullcline 901 u be d c Equilibrium points 00 10 dbc bddcbc c ad c2 Dynamics when a d 1 u nullcline u 0 v 1 uu b be 1 c Equilibrium points 00 10 v nullcline v 0 u be b1 c bc 1 c 1 c2 39 1 1 b and 10 is a sink which attracts any initial values Case 1 cgt Two equilibrium points at 00 and 10 1 Case 2 0 lt c lt 1 H Three equilibrium points at 00 10 uovo bCC7 b1 CC 2bC 1 Stability of 110220 A brief theory of periodic orbits Theorem 1 Inside a periodic orbit there is at least one equilib rium point Corollary 2 Let 0340 be an equilibrium point If there is a stable or unstable orbit of this equilibrium going to infinity as t gt oo then there is no periodic orbit around this equilibrium point Corollary 3 Let 0340 be an equilibrium point If the direction of vector field on the nullclines near 0340 is not clockwise or counterclockwise then there is no periodic orbit around this equilibrium point Theorem 4 Poincar Bendixon If there is a donut region a region without equilibrium point that the vector field is pointing inside at any point on the boundary then there is a periodic orbit in the donut region Corollary 5 If in a region there is only one equilibrium point which is a source or spiral source and the vector field is pointing inside at any point on the boundary then there is a periodic orbit around this equilibrium point A typical case the solutions from outside are spiraling in and the equilibrium point is a spiral source Theorem 6 Poincare Andronov Hopf bifurcation If a bifurca tion occurs such that an equilibrium point changes from a spiral sink to spiral source and the solutions from outside are always spiraling in then there is a periodic orbit around the spiral source Jules Henri Poincare 1854 1912 One of greatest mathematicians in 20th century httpwww gapdcsst andacuk historyMathematiciansPoincareh1 8 du uv do uv 1 d5 M u ub d5 cvub 2 Stability Of u0v0 C1 c bc b C J 1 c 1 c bc O c1 c bc bgt Eigenvalue equation A2 1 c 1 c bc 2 O 1 b S iraI sink or sink 1 b bltO or gt D c c C 1b 1 b Spiral source or source 1 c 96 b gt O or O lt c lt 1 b 9 Conclusion c is where Hopf bifurcation occurs When 1 gt 1b but near 1b is a spiral sink and all solu c u v 1b 1b 0 0 tions tend to 110220 in an oscillating fashion when c passes 1 b but near 1 b 1 b solutions away from 110220 still tends toward 110220 spiral in ward thus a periodic solution with small oscillation emerges around uovo and it is a limit cycle 110220 becomes a spiral source but all Examples b 05 the Hopf bifurcation point is c 13 m 033333 b 02 the Hopf bifurcation point is c 23 m 066667 10 Theorem proved around 1980 For ul u z cv I m d 7 7 s u I b d5 u I b 1 1 b if gt then 10 is as m toticall stable if lt lt 0 1 y I0 y 1 b 6 1 1 1 b then is as m toticall stable and if 0 lt c lt 11 U0 110 y D V 11 then there exists a unique periodic solution and it attracts all positive initial values except the equilibrium points 00 10 and u0 00 Hilbert problem 16 Consider u fuv 12 guv If f and g are polynomials prove the system has only finite many limit cycles When f and g are quadratic functions prove there are at most three limit cycles Counterexample is found by S Shi in 1979 four limit cycles 11 Turing instability i iluum u1 71 Jib 1 6 01 m vvm cvu b 5136 01 4 ux07t u17t vx07t 11361 2 0 uaO u0a Z O vaO 12033 Z O 13 6 01 Nope since 91 O we need fugv lt O for Turing instability 1 b But many Turing type bifurcations occurs when 0 lt c lt 1 b along the equilibrium solution u0v0 The global dynamics of the reaction diffusion system is still not known 12 Neural modeling excitable systems FitzHugh Nagumo system evv al v w dt dw b 5 v c w dt 2225 excitation variable wt potassium conductance a 01 b 05 c 0 e 001 equilibrium 00 is asymptoti cally stable but it is an excitable system where a large excursion occurs if a perturbation is beyond a threshold a 01 b 05 c 01 e 001 equilibrium 010 is unstable and a limit cycle exists relaxation oscillation 13 A review of stability and dynamical behaviors of differential equations scalar ODE at fu system of ODEs at u v vi 9uav7 reaction diffusion equation ut DAu fu 1 E Q with boundary condition reaction diffusion system D A ut u u v 1 E Q with boundary condition vi DyA U gu 2 All equation is in form of Ut FU where u can be a scalar or vector spatial independent or dependent Abstract Equation Ut FU Equilibrium solution U0 such that FUO O linearized operator F UO U0 is stable if the eigenvalues of equation F U0w Aw are all with negative real parts Linear behavior Since the equation Ut FU is approximately the linearized equation Ut F UU U0 near the equilibrium solution U 2 U0 then Ut m ZCiexpOWtwi near U 2 U0 where Moi are the eigenvalue eigenvector pairs of F U0w Aw If A1 max AZ then Ut m 01 expgt1tglgt1 if UO m U0 Example 1 Scalar equation at fu fu0 O linearized operator f u0 f u0 lt O stable otherwise unstable Example 2 Linear system at 2n 31 linearized system J 2 3 vi 4n 31 43 4 3 eigen Dairs A1 1 1 1 1A2 6 2 34 Solution lt gt cletlt 11 gt ch6t lt i ch6t lt 2 if perturbed from equilibrium Eigenvalue problem 2 3 gt10 2 Aw Example 3 Nonlinear ODE system xt2x x2 xy yt 3y 92 2mg Equilibrium points 00 20 30 and 11 2 2x y 13 linearized operator J lt 2y 3 2y 2 For 11 eigenvalues A1 1 1 058 082 A2 1 5 a2 058 082 I 1 058 So near the saddle ponnt 11 lt gt lt 1 016 1 lt O82gt CQe 1 t gt m lt i gt cle 1 tlt 305882 gt when DEF turbed from equilibrium Example 4 scalar diffusion equation ut DUZUZUI a 6077T1t gt 01 ut77r Equilibrium solution u13 O Linearized operator F uw Dwm for w satisfying 1110 2 w7r O Eigenvalue problem Dwm Aw a E 07r 1110 2 w7r O Eigenvalue eigenfunction pairs Am 2 Dm2 ma 2 Sime m 1 2 Solution u13t Zomexp Dm2t Sinmx m c1 exp Dt sinx near u O and u13 O is a stable equilibrium solution Example 5 scalar diffusive Malthusian equation at Dum au 13 E 07r t gt O utO ut7r O Equilibrium solution u13 O Linearized operator F uw Dwm aw for w satisfying 1110 2 w7r O Eigenvalue problem Dwm I aw Aw x E 07r 1110 2 w7r O Eigenvalue eigenfunction pairs Am 2 a D7712 ma 2 Sime m 1 2 Solution u13t Zomexp a Dm2tsinmx m c1 expa Dtsinx near u O u13 O is a stable equilibrium solution if a lt D and it is unstable if a gt D Thus a D is a potential bifurcation point Example 6 Fisher equation at Dum au1 u 13 E 07r t gt O utO ut7r O Equilibrium solution u13 O Linearized operator F uw Dwm aw for w satisfying 100 w7r 0 Eigenvalue problem Dwm I aw Aw x E 07r 1110 2 w7r O Eigenvalue eigenfunction pairs Am 2 a D7712 ma 2 Sime m 1 2 Solution u13t Zomexp a Dm2tsinmx m c1 expa Dtsinx near u O u13 O is a stable equilibrium solution if a lt D and it is unstable if a gt D Thus a D is a bifurcation point Fisher equation cont For a gt D Fisher equation has a positive equilibrium solution uaa from local bifurcation theory and global bifurcation of time mapping uaa is a stable equilibrium solution see lecture notes for a proof Diffusion is a stabilizing process In Fisher equation at Dum au1 u a E 07r t gt O ut O ut7r O Equilibrium solution u O is also an equilibrium solution of at au1 u but an unstable one since f O gt O eigenvalue A elfO gt O Equilibrium solution u13 O is stable for Fisher equation when a lt D and it is still unstable when a gt D eigenvalues Am 2 af 0 Dm2 Thus diffusion makes the eigenvalue more negative and dif fusion in scalar equation has a stabilizing effect If stable for ODE also stable for PDE even unstable for ODE still can be stable for PDE Next model most useful one reaction diffusion system at Duum fuv vi vam guv ut7 ut7 Oi vt7 vt7 Z 0 If u0v0 is an equilibrium solution that fu0vo gu0v0 0 then u0vo is also an equilibrium solution of at fuv vtguv Linearized operator ODE Jacobian J lt f f 9a 911 I duA O fu fv PDE dzagduAadvA J lt 0 dvAv gt lt ya 91 DuAcf fuu07 110M fvu07 110W M Eigenvalue problem DUAgb guu0 UOM gvu0 vow 2 Mb CharaJ 07 O 10
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'