UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A reaction-diffusion model for intracellular calcium Cheung, Robert Wai 1984

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1984_A6_7 C45.pdf [ 2.89MB ]
Metadata
JSON: 831-1.0079345.json
JSON-LD: 831-1.0079345-ld.json
RDF/XML (Pretty): 831-1.0079345-rdf.xml
RDF/JSON: 831-1.0079345-rdf.json
Turtle: 831-1.0079345-turtle.txt
N-Triples: 831-1.0079345-rdf-ntriples.txt
Original Record: 831-1.0079345-source.json
Full Text
831-1.0079345-fulltext.txt
Citation
831-1.0079345.ris

Full Text

A REACTION-DIFFUSION MODEL FOR INTRACELLULAR CALCIUM By ROBERT WAI CHEUNG B . A . S c , The U n i v e r s i t y of B r i t i s h C olumbia, 1982 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE i n THE FACULTY OF GRADUATE STUDIES Department of Mathematics We a c c e p t t h i s t h e s i s as c o n f o r m i n g t o the r e q u i r e d s t a n d a r d THE UNIVERSITY OF BRITISH COLUMBIA October 1984 © Robert Wai Cheung, 1984 86 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the r e q u i r e m e n t s f o r an advance degree a t the U n i v e r s i t y of B r i t i s h C o l umbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g of t h i s t h e s i s f o r s c h o l a r l y purposes may be g r a n t e d by the head of my department or by h i s or her r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . Department of Mat hematics The U n i v e r s i t y of B r i t i s h Columbia 1956 Main M a l l Vancouver, Canada V6T 1Y3 Date October 15, 1984 i i ABSTRACT I n t r a c e l l u l a r c a l c i u m ions ( C a 2 + ) p l a y important r o l e s i n neurobiology by e i t h e r t r i g g e r i n g or modulating a l a r g e number of processes which are a s s o c i a t e d with nerve c e l l behaviour. In s p i t e of t h i s importance, i t i s very d i f f i c u l t e x p e r i m e n t a l l y to o b t a i n q u a n t i t a t i v e i n f o r m a t i o n on the dynamics of i n t r a c e l l u l a r c a l c i u m . For t h i s reason, Connor and Nikolakopoulou [13] formulated a d i f f u s i o n model based upon e x p e r i m e n t a l l y estimated parameters. They used t h i s model to study the s p a t i a l d i s t r i b u t i o n of c a l c i u m w i t h i n the cytoplasm and the i n c r e a s e s i n the c o n c e n t r a t i o n of i n t r a c e l l u l a r c a l c i u m ions that a given i n f l u x of c a l c i u m through the c e l l membrane can b r i n g about. The mathematical model c o n s i s t s of a system of f i v e r e a c t i o n -d i f f u s i o n equations which i s t r e a t e d as a two-point i n i t i a l -boundary-value problem with constant i n i t i a l s t a t e s and n o n l i n e a r boundary c o n d i t i o n s . A n a l y t i c a l r e s u l t s r e v e a l that the system admits a unique s p a t i a l l y homogeneous s t a t i o n a r y s t a t e which i s a s y m p t o t i c a l l y s t a b l e . A r e g u l a r p e r t u r b a t i o n technique i s used f o r c o n s t r u c t i n g an approximate t r a n s i e n t s o l u t i o n which possesses t r a n s i t i o n l a y e r s . A s i m p l i f i e d comparison theorem for p a r a b o l i c equations i s used to p r o v i d e a n a l y t i c a l bounds on the approximate s o l u t i o n . F i n a l l y , a B - s p l i n e c o l l o c a t i o n code, c a l l e d PDECOL, i s used to supplement the a n a l y t i c a l r e s u l t s to provide a d e t a i l e d d e s c r i p t i o n of the c a l c i u m dynamics. i i i TABLE OF CONTENTS pages 1. INTRODUCTION 1 2. THE MATHEMATICAL MODEL 5 3. MATHEMATICAL FORMULATION 13 4. MATHEMATICAL PROPERTIES 17 4.1 EXISTENCE AND UNIQUENESS OF THE STEADY STATE 21 4.2 ASYMPTOTIC STABILITY OF THE STEADY STATE ...28 5. SOME FEATURES OF THE TRANSIENT BEHAVIOUR 34 6. NUMERICAL SOLUTION 48 7. CONCLUSION 56 REFERENCES 58 APPENDIX A 60 APPENDIX B 67 LIST OF TABLES iv pages 1. Parameter values used i n d i f f u s i o n model 12 2. Values of the dimensionless parameters 16 3. Summary of the true maximum values of (U,V,W) together with t h e i r r e s p e c t i v e upper bounds 47 4. Summary of the i n i t i a l , maximum, and s t e a d y - s t a t e c o n c e n t r a t i o n l e v e l s of CaArz, C a 2 + , and CaB 55 V LIST OF FIGURES pages 1. S t r u c t u r e of d i f f u s i o n model with l o c a l r e a c t i o n s . 7 2. Geometrical i l l u s t r a t i o n of the uniqueness of the s t e a d y - s t a t e s o l u t i o n 26 3a-b. Comparison between the exact s o l u t i o n of V with the p e r t u r b a t i o n s o l u t i o n f o r t = 0.001,0.01 ....41 4a. C o n c e n t r a t i o n p r o f i l e s of CaArz d u r i n g a 0.1 second l o a d i n g p e r i o d 49 4b. C o n c e n t r a t i o n p r o f i l e s of C a 2 + d u r i n g a 0.1 second l o a d i n g p e r i o d 50 4c. C o n c e n t r a t i o n p r o f i l e s of CaB d u r i n g a 0.1 second l o a d i n g p e r i o d 51 5a. C o n c e n t r a t i o n p r o f i l e s of CaArz a f t e r i n f l u x 52 5b. C o n c e n t r a t i o n p r o f i l e s of C a 2 + a f t e r i n f l u x 53 5c. C o n c e n t r a t i o n p r o f i l e s of CaB a f t e r i n f l u x 54 6. Sketch of the t r a p e z o i d a l - s h a p e d f u n c t i o n H(t) ....65 CHEMICAL NOTATIONS c o n c e n t r a t i o n of a chemical c a l c i u m i n d i c a t o r dye arsenazo I I I t o t a l c o n c e n t r a t i o n of Arz i n t r i n s i c b u f f e r t o t a l c o n c e n t r a t i o n of B c a l c i u m ions r e s t value of C a 2 + c o n c e n t r a t i o n i n t r a c e l l u l a r c o n c e n t r a t i o n of Ca 2 e x t r a c e l l u l a r c o n c e n t r a t i o n of C a 2 complex formed by C a 2 + and Arz complex formed by C a 2 + and B d i f f u s i o n c o e f f i c i e n t , i = 1,2,3 r e a c t i o n r a t e c o n s t a n t ^ i = 1,2,3, v i i ACKNOWLEDGEMENTS I am indebted to my s u p e r v i s o r , P r o f e s s o r R.M. Miura, f o r suggesting the t h e s i s t o p i c and fo r p r o v i d i n g me with f i n a n c i a l support and p e r c e p t i v e c r i t i c i s m throughout the implementation of t h i s t h e s i s . But most of a l l , I must thank him f o r h i s p a t i e n c e . I must a l s o thank P r o f e s s o r G. Bluman f o r ex p r e s s i n g i n t e r e s t i n my t h e s i s and reading the f i n a l v e r s i o n . This thesis is dedicated to my parents. 1 1. INTRODUCTION In t h i s t h e s i s , the phrase nonlinear diffusion equations r e f e r s s p e c i f i c a l l y t o s e m i l i n e a r systems of second-order p a r t i a l d i f f e r e n t i a l e q u a t i o n s of the form 3u/3t = DAu + F ( x , t , u , V u ) , u = (u ,...,u ). (1.1) 1 m n n n Here x e R , A = E 9 2 / 3 x 2 i s the L a p l a c e o p e r a t o r i n R, and the j = 1 j diffusion matrix D i s d i a g o n a l 1 c o n s i s t i n g of c o n s t a n t n o n n e g a t i v e e l e m e n t s 2 o n l y . D i f f u s i o n i s a p h y s i c a l phenomenon by which m a t t e r i s t r a n s p o r t e d from one p a r t of a system t o ano t h e r as a r e s u l t of random m o l e c u l a r motions. D i f f u s i o n e q u a t i o n s of type (1.1) w i t h a p p r o p r i a t e i n i t i a l and boundary c o n d i t i o n s a r e w i d e l y used as models f o r b i o l o g i c a l , e c o l o g i c a l , and c h e m i c a l systems. In most p h y s i c a l a p p l i c a t i o n s , the independent v a r i a b l e s x and t r e p r e s e n t p o s i t i o n s i n space and t i m e , r e s p e c t i v e l y . In the c o n t e x t of e c o l o g y , the components of u r e p r e s e n t the p o p u l a t i o n d e n s i t i e s of v a r i o u s s p e c i e s , the term DAu r e p r e s e n t s the r a t e of 'In a n i s o t r o p i c media, e.g., c r y s t a l s and t e x t i l e f i b r e s , i n which the m o l e c u l e s have a p r e f e r e n t i a l d i r e c t i o n of o r i e n t a t i o n , D i s g e n e r a l l y not d i a g o n a l . 2 I n many systems, e.g., the i n t e r d i f f u s i o n of me t a l s or. the d i f f u s i o n of o r g a n i c vapours i n high-polymer s u b s t a n c e s , D depends on the c o n c e n t r a t i o n of d i f f u s i n g s u b s t a n c e s , u,,u 2,... In t h i s case and a l s o when the medium i s not homogeneous so t h a t D v a r i e s from p o i n t t o p o i n t , the d i f f u s i o n term i n (1.1) becomes ADu, where D may be a f u n c t i o n of x and u. 2 change of the p o p u l a t i o n d e n s i t i e s at any given p o s i t i o n and time due to random s p a t i a l m i g r a t i o n , and the term F r e p r e s e n t s the i n t e r a c t i o n s between v a r i o u s s p e c i e s such as r e p r o d u c t i o n processes or deaths. V a r i o u s examples of b i o l o g i c a l and e c o l o g i c a l d i f f u s i o n models are given i n [4,12], In t h i s t h e s i s , we study a s p e c i a l case of equation (1.1) which models the dynamics of i n t r a c e l l u l a r c a l c ium ion c o n c e n t r a t i o n changes i n one s p a t i a l dimension. I n t r a c e l l u l a r c a l c i u m e i t h e r t r i g g e r s or modulates a l a r g e number of processes which are a s s o c i a t e d with nerve c e l l behaviour. Some of these are muscle c o n t r a c t i o n , r e l e a s e of n e u r a l t r a n s m i t t e r , and sensory c e l l t r a n s d u c t i o n . The mathematical model [13] c o n s i s t s of a system .of f i v e n o n l i n e a r d i f f u s i o n equations d e s c r i b i n g the e v o l u t i o n of the v a r i o u s i o n i c c o n c e n t r a t i o n s which are taken i n t o account. The system of equations i s of the form 3u/9t = D3 2u/3x 2 + F ( u ) , (1.2) where u = (u, , u 2 , u 3 , u,, u 5) , x e [ a , b ] , D' i s a d i a g o n a l matrix with constant p o s i t i v e e n t r i e s d , and F i s a q u a d r a t i c p o l y -i nomial i n u. The i n i t i a l c o n d i t i o n s are s p e c i f i c c o n s t a n t s chosen to r epresent the e q u i l i b r i u m s t a t e of the system. The boundary c o n d i t i o n s are of n o n l i n e a r Neumann type. Equations of type (1.2) with F polynomials i n u 1,u 2,..,u 5 are sometimes known as react ion-diffusion equations. Although (1.2) i s a s p e c i a l case of (1.1), i t i n c l u d e s two very important extreme cases : 3 i ) F = 0, d > 0 : i 3u /3t = d 3 2u /3x 2, i = 1,2,...,5. (1.3) i i i Thus each equation i s a s c a l a r d i f f u s i o n equation. i i ) d = 0 : i du/dt = F ( u ) , ( 1 .4) which i s a system of k i n e t i c equations a s s o c i a t e d with (1.2). Observe that even when D # 0, equation (1.4) i s s a t i s f i e d by x-independent s o l u t i o n s of (1.2). The layout of t h i s t h e s i s i s as f o l l o w s . In Chapter 2 we in t r o d u c e the mathematical model proposed by Connor and Nikolakopoulou [13] (CN) in t h e i r e f f o r t s to study the dynamics of i n t r a c e l l u l a r c a l c i u m ion c o n c e n t r a t i o n changes i n nerve cytoplasm 1 r e s u l t i n g from e x t e r n a l i n f l u x e s of c a l c i u m i o n s . Numerical s o l u t i o n s to the mathematical model are given i n t h e i r paper. However, they have i n c l u d e d no a n a l y t i c a l r e s u l t s . I t i s our o b j e c t i v e in t h i s t h e s i s to provide some a n a l y t i c a l d e s c r i p t i o n s of the c a l c i u m dynamics and compare them with the numerical s o l u t i o n s . In Chapter 3 we reformulate the CN model so that i t i s amenable to v a r i o u s mathematical analyses and i n Chapter 4 we C y t o p l a s m i s c h a r a c t e r i z e d as a simple d i f f u s i o n compartment c o n t a i n i n g an i n t r i n s i c b u f f e r f o r Ca. 4 explo r e some of i t s mathematical p r o p e r t i e s . These p r o p e r t i e s are used i n subsequent s e c t i o n s to f a c i l i t a t e the c a l c u l a t i o n s of the s t e a d y - s t a t e s o l u t i o n and to e s t a b l i s h i t s uniqueness and asymptotic s t a b i l i t y . In Chapter 5, we use a r e g u l a r p e r t u r b a t i o n technique to obt a i n an approximate small-time s o l u t i o n to the model equations so that we can make a q u a n t i t a t i v e statement about the t r a n s i e n t behaviour of the cal c i u m c o n c e n t r a t i o n . Moreover by u t i l i s i n g a comparison theorem, we provide a n a l y t i c a l bounds on the exact s o l u t i o n . In Chapter 6 we present the numerical s o l u t i o n s to the problem and compare them with the experimental r e s u l t s given by Connor and Nikolakopoulou [13]. In Appendix A, we give a b r i e f d e s c r i p t i o n of the numerical a l g o r i t h m , PDECOL, which was used to so l v e the model equations and r e p o r t some of the d i f f i c u l t i e s encountered i n the numerical procedures. In Appendix B, we f i n d an approximate s o l u t i o n to a l i n e a r second-order p a r t i a l d i f f e r e n t i a l equation with n o n l i n e a r boundary c o n d i t i o n s and the r e s u l t s are used to s i m p l i f y the c a l c u l a t i o n s f o r the t r a n s i e n t s o l u t i o n s . A l l computations l e a d i n g to the numerical r e s u l t s presented in t h i s paper were performed i n double p r e c i s i o n at the U n i v e r s i t y of B r i t i s h Columbia Computing Center using an Amdahl V/8 Model 470 computer. 5 2. THE MATHEMATICAL MODEL There i s a c o n s i d e r a b l e amount of evidence i n the l i t e r a t u r e (Sandow [23], J o b s i s & O'Connor [18], Ebashi & Endo [16]) to suggest that small t r a n s i e n t changes i n sarcoplasmic c a l c i u m c o n c e n t r a t i o n play a c r u c i a l p a r t i n the complex process of e x c i t a t i o n - c o n t r a c t i o n c o u p l i n g i n s k e l e t a l muscle. I t i s thus of i n t e r e s t to know the p r e c i s e s i z e s of these c a l c i u m changes that l e a d to the pro d u c t i o n of t e n s i o n . Various experimental methods such as photometric [ 6 ] , v o l t a g e clamp technique [ 1 ] , and i n d i c a t o r dye absorbance [9] have been designed f o r t h i s p a r t i c u l a r purpose. Among them, i n d i c a t o r methods have proven t o be the most s u c c e s s f u l because of t h e i r s e n s i t i v i t y and a b i l i t y to respond r a p i d l y to c o n c e n t r a t i o n changes. However, such methods are imprecise as to the l o c a t i o n of changes w i t h i n the cytoplasm. S e v e r a l reviews of the i n d i c a t o r techniques are given i n [8] and [ 9 ] . Experimental r e s u l t s by Ahmed and Connor [2] using i n d i c a t o r techniques show that the i n t e r n a l c a l c i u m c o n c e n t r a t i o n ( [ C a 2 + ] ) changes are brought about e i t h e r by a c a l c i u m i n f l u x from some e x t r a c e l l u l a r r e s e r v o i r or by the r e l e a s e of ca l c i u m ions from an i n t e r n a l source. However, the short time course of the ca l c i u m i n f l u x (on the order of ten m i l l i s e c o n d s ) and the small c o n c e n t r a t i o n changes of C a 2 + ( i n the micromolar range) together with the small volumes g e n e r a l l y encountered with b i o l o g i c a l p r e p a r a t i o n s make the problem of d i r e c t q u a n t i t a t i v e measurement of the C a 2 + t r a n s i e n t a very d i f f i c u l t one. 6 In 1982 Connor and Nikolakopoulou performed experiments 1 on molluscan g i a n t neurons using i n d i c a t o r dye absorbance in combination with a v o l t a g e clamp technique i n an attempt to make a q u a n t i t a t i v e statement about the dynamics of i n t r a c e l l u l a r c a l c i u m ion c o n c e n t r a t i o n changes. The experiments p r o v i d e d them with a great deal of q u a l i t a t i v e i n f o r m a t i o n about c a l c i u m e n t r y i n t o the neurons but with l i t t l e or no i n s i g h t as to the s p a t i a l d i s t r i b u t i o n of the C a 2 + w i t h i n the cytoplasm. To add to t h e i r f r u s t r a t i o n , the experiment provided no i n f o r m a t i o n as to the magnitude of the p e r t u r b a t i o n i n [ C a 2 + ] that a given i n f l u x can b r i n g about. In order to s o l v e the p u z z l e , Connor and Nikolakopoulou formulated a o n e - d i m e n s i o n a l 2 d i f f u s i o n model based upon e x p e r i m e n t a l l y e s t i m a t e d 3 parameters and what appear to be reasonable estimates f o r those parameters which they have not been able to measure. F i g u r e 1 i n d i c a t e s the s t r u c t u r e of the CN d i f f u s i o n model and shows that e x t r a c e l l u l a r c a l c i u m ions, C a 0 , are loaded through the membrane f o r a p e r i o d t * lOOmsecs). As soon as the c a l c i u m ions enter the c e l l membrane, they are assumed to be i n s t a n t a n e o u s l y at l o c a l e q u i l i b r i u m with the i n t r i n s i c ' D e t a i l e d d e s c r i p t i o n s of the experimental methods are given in Ahmed & Connor [ 1 ] , Ahmed, et a l . [ 2 ] . 2 I t w i l l soon be c l e a r from the model equations that a more r e a l i s t i c model in 3-dimensions would be mathematically too complicated to handle. 3 E x p e r i m e n t a l methods used in e s t i m a t i n g the parameters of the model are d e s c r i b e d i n [13]. 7 b u f f e r ( B ) 1 , and the c a l c i u m i n d i c a t o r dye arsenazo III ( A r z ) 2 . I n i t i a l l y , the two b u f f e r s (B and Arz) have uniform c o n c e n t r a t i o n s and thus the e n t r y of c a l c i u m ions causes d i f f u s i o n a l movement w i t h i n the membrane. x=0 e x t r a c e l l u l a r space C a 0 s u r f a c e membrane i n t r a c e l l u l a r space v I r 1 1 CaArz Arz + Ca 2*| + B CaB . I < ' , < 1 I ki L i k„ I i i L i T r J V V V V v D, D, D 3 D 2 D 2 D i f f u s i o n a l movement along the x - d i r e c t i o n F i g u r e 1 : S t r u c t u r e of d i f f u s i o n model with l o c a l r e a c t i o n s . Let U , ( x , t ) , U 2 ( x , t ) , U 3 ( x , t ) , U„(x,t), and U 5 ( x , t ) 1An i n t r i n s i c b u f f e r i s an inhe r e n t substance i n the s o l u t i o n that makes the degree of a c i d i t y (hydrogen ion c o n c e n t r a t i o n ) r e s i s t a n t to change when an a c i d or base i s added. 2The absorbance spectrum of the dye arsenazo I I I undergoes known changes when the dye forms a complex with c a l c i u m . These changes have been used i n a number of s t u d i e s to d e s c r i b e c a l c i u m movement l e v e l s i n l i v i n g c e l l s . (See Baker, Hodgkin & Ridgway [5].) 8 represent the c o n c e n t r a t i o n s of CaArz, Arz, C a 2 + , B, and CaB, r e s p e c t i v e l y . R e c a l l from elementary chemistry that the law of mass action f o r a f i r s t - o r d e r chemical r e a c t i o n between species A, B, and C, k _7 B + C i s given by d[A]/dt = k-[B][C] - k + [ A ] , d[B]/dt = k + [ A ] - k- [ B ] [ C ] , d[C]/dt = k + [ A ] - k- [ B ] [ C ] , where [•] denotes the c o n c e n t r a t i o n of a chemical and where k + and k" denote the forward and reverse r e a c t i o n r a t e constants of chemical A, r e s p e c t i v e l y . Thus assuming a f i r s t - o r d e r r e a c t i o n i n an i s o t r o p i c medium and i g n o r i n g s p a t i a l i n h o m o g e n i t i e s 1 , we can apply the law of mass a c t i o n and write down the k i n e t i c equations plus the d i f f u s i o n terms a s s o c i a t e d with our chemical model : 3U,/9t = D , 3 2U 1 / 3 x 2 + k,U 2U 3 - k 2U,, (2 .1a) 3U 2 /3t = D 1 3 2 U 2 / 3 x 2 - k!U 2U 3 + k 2U,, (2.1b) 3U 3 /3t = D 3 3 2 U 3 / 3 x 2 - k,U 2U 3 + k 2U, - k 3 U 3 U „ + k„U 5, (2 .1c) 3U„/3t = D 23 2U«/3x 2 - k 3 U 3 U „ + k u U 5 , (2 .Id) 3U 5 /3t = D 2 3 2 U 5 / 9 x 2 + k 3 U 3 U 4 - k f lU 5, (2 .1e) where D 1 f D 2, and D 3 are constant d i f f u s i o n c o e f f i c i e n t s , and k 1 f 'These are reasonable assumptions to ensure that the d i f f u s i o n matrix D i s di a g o n a l and c o n s t a n t . 9 k 2, k 3, and k„ are the forward and reverse r e a c t i o n r a t e c o n s t a n t s as shown i n F i g u r e 1. It i s c l e a r from the chemistry of our p h y s i c a l model that [Arz] and [B] i n the system are conserved at each x and f o r a l l t > 0; we s h a l l denote these c o n c e n t r a t i o n s by [Arz] and [B] , T T r e s p e c t i v e l y . But moreover, s i n c e the i n i t i a l d i s t r i b u t i o n i s uniform and s i n c e U, and U 2 have the same d i f f u s i o n c o e f f i c i e n t , we have [Arz] = U,(x,t) + U 2 ( x , t ) , (2.2a) T fo r t > 0 and a l l x. A l t e r n a t i v e l y , (2.2a) can be d e r i v e d by adding equations (2.1a) and (2.1b) g i v i n g 3(U,+U 2)/3t = D 1 3 2 ( U 1 + U 2 ) / 3 x 2 . Now s i n c e there i s no f l u x c o n t r i b u t i o n from U, and U 2 at the boundaries and s i n c e (2.2a) holds at t = 0, i t holds f o r a l l t > 0. S i m i l a r l y [B] = U«(x,t) + U 5 ( x , t ) . (2.2b) T Thus equations (2.2a)-(2.2b) allow us to e l i m i n a t e two of the f i v e equations i n (2.1a-e). We choose to e l i m i n a t e U 2 and U*, s i n c e they correspond to the c o n c e n t r a t i o n s of the b u f f e r s and hence are of l e s s i n t e r e s t . The reduced system becomes 3U,/3t = D 13 2U,/3x 2 + k, ( [ A r z ] - U,)U 3 - k 2 U 1 f T 3U 3/3t = D 33 2U 3/3x 2 - k , ( [ A r z ] - U,)U 3 + k 2U! T - k 3 ( [ B ] - U 5 ) U 3 + k a U 5 , T 3U 5/3t = D 23 2U 5/3x 2 + k 3 ( [ B ] - U 5 ) U 3 - k f tU 5. (2.3) T 10 In order to completely s p e c i f y the problem, we must impose a p p r o p r i a t e i n i t i a l and boundary c o n d i t i o n s to (2.3). To t h i s end observe that the system i s i n i t i a l l y at e q u i l i b r i u m and thus the i n i t i a l c o n d i t i o n s can be obtained by simply s e t t i n g the k i n e t i c terms to zero, g i v i n g U,(x,0) = k , U 3 ( x , 0 ) [ A r z ] / ( k 2 + k,U 3(x,0)) i M / 1 , T U 3(x,0) = [ C a 2 + ] uM/i R U 5(x,0) = k 3 U 3 ( x , 0 ) [ B ] /(k, + k 3 U 3 ( x , 0 ) ) M M / 1 , (2.4) T where [ C a 2 + ] denotes the i n i t i a l r e s t value of c a l c i u m ion R c o n c e n t r a t i o n and yM/1 stands f o r micromolars per l i t e r . The boundary c o n d i t i o n s must take i n t o account the i n f l u x of C a 2 + at x = 0. By F i c k ' s law [15], we have 9U 3(0,t)/9x = I / ( z F D 3 ) , Ca where I = c a l c i u m c u r r e n t , z = valence of C a 2 + = 2, Ca F = faraday constant (96500 coulombs/mole). From Ohm's Law I = g V, Ca Ca where g = conductance of membrane, Ca V = net d r i v i n g p o t e n t i a l on C a 2 + = E - E , m Ca where E i s the membrane p o t e n t i a l and E i s the Nernst m Ca 11 p o t e n t i a l f o r C a 2 + [21] given by E = ( 2 . 3 0 3 R T / z F ) l o g 1 o [ C a 2 + ] / [ C a 2 * ] V o l t s , Ca 0 i where R = u n i v e r s a l gas constant (8.31 joules/mole/°K), T = temperature in °K, [ C a 2 + ] = the e x t r a c e l l u l a r c o n c e n t r a t i o n of C a 2 + , 0 [ C a 2 + ] = the i n t r a c e l l u l a r c o n c e n t r a t i o n of C a 2 + . i Thus assuming a temperature of 18°C, the Nernst p o t e n t i a l f o r C a 2 + i s E = 2 9 l o g 1 0 [ C a 2 + ] / [ C a 2 + ] mVolts. Ca 0 i F i n a l l y , s i n c e a l l other f l u x e s are zero at the boundaries, the a p p r o p r i a t e boundary c o n d i t i o n s f o r our problem are 3U,(0,t)/3x = 0 = 3U,(L,t)/3x, 3U 5(0,t)/3x = 0 = 3U 5(L,t)/3x, 3U 3(0,t)/3x = (g / z F D 3 ) ( E - 2 9 1 o g , 0 [ C a 2 + ] 0 / U 3 ( 0 , t ) ) yM/lm, Ca m 3U 3(L,t)/3x = 0, where L r e p r e s e n t s the r a d i u s of a neuron. For the purpose of mathematical c a l c u l a t i o n s , we can choose L s u f f i c i e n t l y l a r g e so that d u r i n g the e n t i r e l o a d i n g p e r i o d t * , the i n f l u x at the boundary (x = 0) does not induce c o n c e n t r a t i o n changes near the center of the neuron. In Table 1, we summarize the parameter values used by Connor and Nikolakopoulou's [13] i n t h e i r d i f f u s i o n model. Numerical r e s u l t s p u b l i s h e d i n Connor and Nikolakopoulou's paper [13] i n d i c a t e that the proposed model giv e s a good account to the t r u e q u a l i t a t i v e behaviour of c a l c i u m dynamics. However, 1 2 the paper c o n t a i n s no mathematical a n a l y s i s of the model equations. I t i s our i n t e n t i o n to adopt the model and use i t to giv e some a n a l y t i c a l d e s c r i p t i o n s of ca l c i u m dynamics. Table 1 : Parameter Values used i n D i f f u s i o n Model (Taken from Connor and Nikolakopoulou [13]) Cytoplasmic and I n d i c a t o r Parameters : [Arz] = 300 M M / 1 , T D, = 0.06 Mm 2/msec, [B] =200 M M / 1 , T D 2 = 0.01 Mm 2/msec, D 3 = 0.2 Aim 2/msec, k, = 0.025 1/MMmsec, k 2 = 1.0/msec, k 3 = 1.0 1/MMmsec, k„ = 1.0/msec. Parameter Values used f o r Calcium I n f l u x : [Ca] = 10000 M M / 1 , 0 [Ca] = 0.1 M M / 1 , R F = 96500 coulombs /M , g = 8.6x10-" mhos/cm2, Ca E = 10 mVolts, m z = 2. C o n t r o l Parameters : t * = 100 msec, L = 30 Mm. Dimensions of the V a r i a b l e s [ x] = am, [t ] = msec, [U, ] = [U 2] = [U 3] = [u 4 ] = t u 5 ] = MM/1. 1 3 3. MATHEMATICAL FORMULATION Exact s o l u t i o n s to d i f f e r e n t i a l equations are r a r e i n p r a c t i c e because of n o n l i n e a r i t i e s , inhomogeneities, and gen e r a l boundary c o n d i t i o n s . As a r e s u l t , s o l u t i o n s are approximated using a n a l y t i c a l techniques, numerical techniques, or combinations of both. In t h i s and the next two cha p t e r s , we s h a l l c o n c e n t r a t e on a n a l y t i c a l techniques, which when combined with the numerical s o l u t i o n d i s c u s s e d i n Chapter 6, give a thorough d e s c r i p t i o n of the behaviour of the model equations (2.3)-(2.5). Our f i r s t step i n s o l v i n g the present problem i s to nondimensionalize the governing equations (2.3)-(2.5); f o r i t i s only i n di m e n s i o n l e s s form that we can compare the order of magnitude of the d i f f e r e n t v a r i a b l e s and parameters of the system and thus be ab l e to i d e n t i f y the important ones. For our problem the d i f f u s i o n c o e f f i c i e n t s s a t i s f y D 3 >> D 2 > D, and thus U 3 d i f f u s e s f a s t e r than U, and U 5. Moreover s i n c e the system i s i n i t i a l l y at e q u i l i b r i u m ( i . e . , the forward and backward r e a c t i o n s are i n b a l a n c e ) , the d i f f u s i o n process of U 3 w i l l , at l e a s t i n e a r l y stages of development, dominate the e n t i r e r e a c t i o n - d i f f u s i o n p r o c e s s . Thus we r e s c a l e the governing equation of U 3 so that the 9/9t and 3 2/9x 2 terms are of the same or d e r . In view of (2.2a)-(2.2b), the other two dependent v a r i a b l e s U, and U 5 can be made dimensionless by n o r m a l i z i n g with re s p e c t to [B] and [Arz] , r e s p e c t i v e l y . The independent time T T v a r i a b l e t can be made dimensionless by using i/k2k„ which has the convenient value of 1. The dimensionless independent v a r i a b l e 1 4 corresponding to x i s now determined. In summary we have introduced the f o l l o w i n g dimensionless independent v a r i a b l e s , x' and t ' , and the dimensionless dependent v a r i a b l e s U, V and W : where t' = / ( k 2 k , ) t , U = U,/[Arz] , T W = U 5/[B] , T D = v/D 3/i/(k 2k 4 ) , x' = x/D, V = kU 3, k = v / ( k , k 3 / k 2 k a ) . (3. 1a) Dimensionless parameters are a , = i/(k ,/k 3 ) , D 1 3 = D,/D3, 6 = k 2 / ( a 2 k , [ A r z ] ), T a 2 = |/( k 2/k„) , D 2 3 = D 2/D 3, e = k„ a 2 / ( k 3 [ B ] ). (3.1b) T A f t e r dropping the primes, equations (2.3)-(2.5) become 3U/3t = D 1 3 3 2 U / 3 x 2 + f( U , V ) , 3V/3t = 3 2V/3x 2 + g(U,V,W), 3W/3t = D 2 33 2W/3x 2 + h(V,W), (3.2) where f(U,V) = a,(1-U)V - a 2U, h(V,W) = (1-W)V/a, - W/a2, g(U,V,W) = - [ ( 1 / ( 6 a , ) ) f ( U , V ) + (a,/t)h(V,W)], subject to the boundary c o n d i t i o n s f o r t > 0, 3U(0,t)/3x = 3W(0,t)/3x = 0, 3U(L,t)/3x = 3V(L,t)/3x = 9W(L,t)/3x = 0, 15 r A + B l n V ( 0 , t ) , 0 < t < t * , 3V(0,t)/3x = 0, t > t * , (3.3) where A = £[E - ( 2 9 / l n l O ) l n ( k [ C a ] 0 ) ] , m B = 29$/ln1(D, $ = (10 7g/zFD 3)kD, and i n i t i a l c o n d i t i o n s on 0 < x < L, U(x,0) = U 0 = V 0 / ( a 2 / a , + V 0 ) , V(x,0) = V 0 = 0.1k, W(x,0) = W0 = V 0 / ( a 1 / a 2 + V 0 ) . (3.4) Hence the present problem depends on s i x parameters, a,, a 2 , e, 6, D 1 3, and D 2 3 . For the set of parameter values given i n Table 1, 0 < e << 5 << 1. Thus s i n c e e and 6 are small we can use them as p e r t u r b a t i o n parameters to o b t a i n an approximate s o l u t i o n to the problem (see Chapter 5). In the s p e c i a l case when e = 0, the s o l u t i o n to (3.2) s u b j e c t to i n i t i a l c o n d i t i o n (3.4) i s simply Note that i f V 0 i s chosen such that l n V 0 = -A/B, the boundary c o n d i t i o n s reduce to homogeneous Neumann types at t = 0. In other words, the r e a c t i o n s are i n a steady s t a t e and the c o n c e n t r a t i o n l e v e l of each chemical w i l l remain at i t s r e s t v a l u e . Henceforth we assume that l n V 0 # -A/B and work with the d i m e n s i o n l e s s equations (3.2)-(3.4) which we s h a l l r e f e r to as system (P). Using the set of parameter values given in Table 1, U(x,t) = U V(x,t) = V 0, W(x,t) = W o • 1 6 we can c a l c u l a t e the dimensionless parameter values from equations (3.1a-b), (3.3) and (3.4). The r e s u l t s are summarized in Table 2. Table 2 : Values of the dimensionless parameters a, = 0.1581 , a 2 = 1.0000, = 0.300, D 2 3 = 0.050, € = 0.0050, 6 = 0 .1333 , D = 0.4472, k = 0.1581 , A = -1.30398, B = 0.19842, . U 0 = 0.00249, V 0 = 0.01581 , W0 = 0.09091. In Chapter 5 we s h a l l r e s c a l e the i n i t i a l v a l u e s U 0, V 0, and W0 to 0(1) before using a r e g u l a r p e r t u r b a t i o n technique to ob t a i n an approximation to the i n i t i a l t r a n s i e n t s o l u t i o n . We s h a l l , however, f i r s t i n v e s t i g a t e some mathematical p r o p e r t i e s of system ( P). 1 7 4. MATHEMATICAL PROPERTIES In the a n a l y s i s of a given problem, i t i s always u s e f u l to f i r s t e x p l o i t the s p e c i a l s t r u c t u r e s of that problem, e.g., symmetry. These s p e c i a l s t r u c t u r e s o f t e n f a c i l i t a t e c a l c u l a t i o n s and, more imp o r t a n t l y , p r o v i d e g r e a t e r i n t u i t i o n and understanding of the problem. In t h i s Chapter, we s h a l l examine some mathematical p r o p e r t i e s of system (P) [equations ( 3 . 2 ) - ( 3 . 4 ) ] . These p r o p e r t i e s are u s e f u l f o r our analyses i n l a t e r c h a p t e r s . Property 1 : S p e c i a l f e a t u r e s of the k i n e t i c terms ' The v e c t o r - v a l u e d f u n c t i o n F = (f,g,h) possesses two s p e c i a l f e a t u r e s which w i l l g r e a t l y s i m p l i f y our l a t e r a n a l y s e s . Before s t a t i n g them, we need the f o l l o w i n g d e f i n i t i o n s : 1. A v e c t o r - v a l u e d f u n c t i o n S(u) i s s a i d to be quas i-monol one nondecreasi ng i n u i f each G i i s nondecreasing i n u , f o r j # i . j 2. For any given V > 0, d e f i n e the set T(V) = {(U,V,W) : U = a , V / [ a 2 + a,V], W = a 2 V / [ a , + a 2 V ] } . The two s p e c i a l f e a t u r e s of F are : 51 : F i s quasi-monotone nondecreasing i n u = (U,V,W). 52 : F(U,V,W) = (f(U,V),g(U,V,W),h(V,W)) = U , f o r a l l V > 0 and (U,V,W) e T ( V ) . S1 can be v e r i f i e d d i r e c t l y by d i f f e r e n t i a t i n g F with 18 respect to u and n o t i n g that both U and W must l i e between 0 and 1 ( c f . Property 3). To prove S2, we simply set f(U,V) = 0 = h(V,W). S2 s t a t e s that our chemical system has a continuum of e q u i l i b r i u m s t a t e s and i n the s p e c i a l case when V = V 0, i t reduces to the statement that the chemical system i s i n i t i a l l y at e q u i l i b r i u m . Property 2 : Conservation Law As i n the case of many p h y s i c a l systems, the s o l u t i o n to system (P) i s governed by a c o n s e r v a t i o n law. To show t h i s observe from equation (3.2) that (l/5a,)3U/3t + 3V/3t + (a,/e)3W/3t = (D, 3/5a,)3 2U/3x 2 + 3 2V/3x 2 + (D 2 3a,/e)3 2W/3x 2. (4.1) I n t e g r a t i n g both s i d e s of equation (4.1) with respect to x from 0 to L and u s i n g boundary c o n d i t i o n s (3.3) we have L J 3z/3t dx = -3V(0,t)/3x, (4.2) 0 where z ( x , t ) = 11/(63,) + V + (a,/e)W. I n t e g r a t i n g equation (4.2) with respect to t , we have L r 0 ( t ) , t < t * , J [ z ( x , t ) - z ( x , 0 ) ] dx = 0 (4.3) P i t * ) , t > t * . where t j3(t) = -/ 3V(0,s)/3x ds. 0 (4.4) 19 Equation (4.3) i s the c o n s e r v a t i o n law governing system (P). Observe that s i n c e [ C a 2 + ] i s loaded i n t o the system, 3 V / ( 0 , t ) / 3 x < 0 for 0 < t < t * , and hence 0 ( t ) i s p o s i t i v e and s t r i c t l y i n c r e a s i n g for 0 < t < t * . In S e c t i o n 4.1 we s h a l l show that the s t e a d y - s t a t e s o l u t i o n to system (P) i s s p a t i a l l y homogeneous and i n such case (4.3) s i m p l i f i e s to z(x,») = z ( x , 0 ) + <3(t*)/L. (4.5) Equation (4.5) s t a t e s that the t o t a l e q u i l i b r i u m c o n c e n t r a t i o n of the chemicals i s the sum of the t o t a l i n i t i a l c o n c e n t r a t i o n plus the t o t a l amount of [ C a 2 + ] i n j e c t e d i n t o the system. We s h a l l see in subsequent s e c t i o n s that equations (4.3)-(4.5) p l a y key r o l e s i n e s t a b l i s h i n g the uniqueness and s t a b i l i t y of the e q u i l i b r i u m s o l u t i o n . Property 3 : I n v a r i a n t Region From the chemistry of the system, one would c o n j e c t u r e that U, V, and W are bounded f u n c t i o n s of x and t . In p a r t i c u l a r , i t i s c l e a r from (2.2a-b) and (3.1a) that U and W must l i e between 0 and 1. We s h a l l now g i v e a mathematical statement of our c o n j e c t u r e which r e l i e s on the concept of an invariant region o r i g i n a l l y i n troduced by Chueh, Conley, and Smoller [ 1 1 ] , An i n v a r i a n t region Z i s a subset of the phase space ( i . e . (U,V,W) space) such that i f the v a l u e s of the s o l u t i o n of the PDE are c o n t a i n e d in the r e g i o n f o r some value of t , e.g., the i n i t i a l c o n d i t i o n s , then the s o l u t i o n remains i n the region f o r 20 a l l l a t e r times. We s h a l l show that an i n v a r i a n t region to system (P) i s given by L = { (U,V,W) : 0 < U 0 ^ U < U < 1 , 0 < V o < V < V <oo, m m 0 < W 0 < W < W < 1 } , (4.6) m where (U ,V ,W ) e T(V ) and V i s a f i n i t e upper bound of V to m m m m m be determined. Since the second s p a t i a l d e r i v a t i v e i s always nonnegative ( n o n p o s i t i v e ) at a minimum (maximum), by the c l o s i n g remark i n Se c t i o n 2 of [14], i t s u f f i c e s to show that F does not po i n t out of I, i . e . , (1) f ( U 0 , V ) > 0, f(U ,V) < 0; m (2) g(U,V0,W) > 0, g(U,V ,W) < 0, m (3) h(V,W 0) > 0, h(V,W ) < 0; m hold f o r a l l (U,V,W) e I . To t h i s end observe that f ( U 0 , V ) = f ( U 0 , V 0 ) + a,(1-U 0 ) ( V - V 0 ) = a,(1-U 0)(V-V 0) ;> 0, and f ( U ,V) = f(U ,V ) + a , ( l - D )(V-V ) = a,(1-U )(V-V ) < 0, m m m m m m m sinc e by S2 f ( U 0 , V 0 ) = 0 = f ( U ,V ). The proof f o r (3) i s m m s i m i l a r . To prove (2) simply observe that g(U,V0,W) > g(U 0,V 0,W 0) = 0, and g(U,V ,W) < g(U ,V ,W ) = 0, m m m m where we have again a p p l i e d S2 i n o b t a i n i n g the right-hand e q u a l i t y . In Chapter 5 we s h a l l o b t a i n an a n a l y t i c a l estimate for V . m 21 4.1 EXISTENCE AND UNIQUENESS OF THE STEADY STATE In t h i s t h e s i s , a ve c t o r u i s d e f i n e d as a steady-state s o l u t i o n to a system of PDEs i f f u i s a time independent s o l u t i o n to the PDEs s a t i s f y i n g the boundary c o n d i t i o n s 1 . For a general system of r e a c t i o n - d i f f u s i o n equations, two que s t i o n s concerning the s t e a d y - s t a t e of the system a r i s e n a t u r a l l y . F i r s t , does there e x i s t a unique s t e a d y - s t a t e s o l u t i o n ? Second, which s t e a d y - s t a t e s o l u t i o n s are s t a b l e ? In t h i s s e c t i o n we s h a l l answer the f i r s t q u e s t i o n and show that f o r a l l p o s i t i v e parameter values a,, a 2 , 5 and e, system (P) has a unique s t e a d y - s t a t e v a l u e . In the next s e c t i o n we show that the s t e a d y - s t a t e s o l u t i o n i s a s y m p t o t i c a l l y s t a b l e . Before p r o v i n g the main r e s u l t , we s h a l l use the concept of i n v a r i a n t r e g i o n developed i n Chapter 4 (Property 3) and s t a t e a lemma which i s a d i r e c t consequence of a theorem proved by Conway, Hoff, and Smoller (Theorem 3.1 i n [1 4 ] ) . Lemma Let u = (u ,...,u ), v = (v ,...,v ), m > 1, x e I = [a , b ] , 1 m 1 m and D be a di a g o n a l matrix with p o s i t i v e e n t r i e s d . Let (S) and i ( S 1 ) denote the problems ^ e have excluded the steady o s c i l l a t o r y s o l u t i o n i n our d e f i n i t i o n . 2 2 Ts) (S') 3u/3t = D 3 2 u / 3 x 2 + F(u) u ( x , t 0 ) = u 0 ( x ) , x e I , 3 u ( x , t ) / 3 x = 0 , x e 31, t * t 0 dv/dt = F(v), b v ( t 0 ) = [/ u 0 ( x ) d x ] / ( b - a ) . F i n a l l y l e t o denote t h e q u a n t i t y a = dX - M where d = min{d }, i i X = s m a l l e s t p o s i t i v e e i g e n v a l u e of - d 2 / d x 2 on I s u b j e c t t o homogeneous Neumann boundary c o n d i t i o n s = j r / ( b - a ) , M = max {|dF(u)| : u e Z}, where I i s the i n v a r i a n t r e g i o n d e f i n e d by e q u a t i o n ( 4 . 6 ) , and the q u a n t i t y |dF(u)| denotes t h e d e t e r m i n a n t of the J a c o b i a n m a t r i x w i t h e n t r i e s 3F /3u . Now i f u 0 ( x ) e £ f o r a l l x c I , and i j a > 0, then u ( x , t ) -> v ( t ) u n i f o r m l y and e x p o n e n t i a l l y f o r a l l x e I as t —> ». The above lemma s t a t e s t h a t t h e s o l u t i o n t o problem (S) can be a p p r o x i m a t e d a r b i t r a r i l y c l o s e l y by t h e s o l u t i o n t o the r e l a t e d k i n e t i c system ( S 1 ) i f t i s s u f f i c i e n t l y l a r g e . I n o t h e r words, when a > 0, the l a r g e time b e h a v i o u r of the s o l u t i o n t o system (S) i s d e t e r m i n e d o n l y by the r e a c t i o n mechanism F ( u ) . 23 In view of the above lemma, we make the f o l l o w i n g p r o p o s i t i o n : P r o p o s i t i o n I : A l l s o l u t i o n s to system (P) decay e x p o n e n t i a l l y to a constant s t e a d y - s t a t e s o l u t i o n f o r any (U 0,V 0,W 0) e I . Proof : For t > t * , the n o n l i n e a r boundary c o n d i t i o n to system (P) reduces to the homogeneous Neumann c o n d i t i o n . By assumption, (U 0,V 0,W 0) e Z and i t f o l l o w s from Property 3 that ( U (x,t*),V(x,t*),W(x,t*)) e Z. Hence the q u a n t i t i e s L u 0 ( t * ) = / U(x,t*)dx/L, 0 L v 0 ( t * ) = / V(x,t*)dx/L, 0 L w 0(t*) = / W(x,t*)dx/L, 0 a l s o belong to Z s i n c e they merely correspond to the s p a t i a l averages of U ( x , t * ) , V ( x , t * ) , and W(x,t*), r e s p e c t i v e l y . Thus i f we choose u 0 ( t * ) , v 0 ( t * ) , and w 0(t*) as the i n i t i a l c o n d i t i o n s to the r e l a t e d k i n e t i c system for t > t * , then by our lemma, we need only to show that a > 0 and that the s t e a d y - s t a t e s o l u t i o n to the r e l a t e d k i n e t i c system i s a c o n s t a n t . We c a l c u l a t e 24 |dF| = -a a/(6a,) 0 a, b -b/5 - c/e c/a , (a,/e)d -d where a = a,V + a 2 , b = (1-U), c = (1-W), d = V/a, + l / a 2 . Observe that the second row i s a l i n e a r combination of the f i r s t and t h i r d rows and thus |dP| = 0 f o r a l l (U,V,W) and a l l parameter v a l u e s a 1 f a 2 , 8, and e. Hence 1 a = dX > 0. Next we s h a l l show that the s t e a d y - s t a t e s o l u t i o n to the r e l a t e d k i n e t i c system i s a cons t a n t . Let z ( t ) = (U-U ,V-V ,W-W ) T c c c denote the d e v i a t i o n from the c r i t i c a l p o i n t (U ,V ,W ) of the c c c k i n e t i c system. By Theorem (9.2) of [ 7 ] , i t s u f f i c e s to c o n s i d e r only the l i n e a r i z e d form of the k i n e t i c system about (U ,V ,w ), c c c z ' ( t ) = A 0 z ( t ) , and show that a l l eigenvalues of the l i n e a r i z e d matrix A 0 are e i t h e r zero or have negative r e a l p a r t s . C l e a r l y A 0 i s the same 1One can i n t e r p r e t the c o n d i t i o n that a be p o s i t i v e to mean that the minimum r a t e of d i f f u s i o n i n the s p a t i a l domain i s la r g e r e l a t i v e to the maximum ra t e of r e a c t i o n . 25 as the Jacobian matrix dF with (U,V,W) re p l a c e d by (U ,V ,W ) i n c c c the d e f i n i t i o n of a, b, c, and d. We s h a l l show in S e c t i o n 4.2 (Property 4, with n = 0) that one eigenvalue of A 0 i s zero and the other two eigenvalues have negative r e a l p a r t s . Thus the s t e a d y - s t a t e s o l u t i o n to the r e l a t e d k i n e t i c system i s a constant and t h i s completes the p r o o f . We now turn to f i n d i n g the s t e a d y - s t a t e s o l u t i o n (U ,V ,W ) s s s of system ( P ) . By P r o p o s i t i o n I, (U ,V ,W ) are independent of x, s s s hence they must s a t i s f y f ( U ,V ) = 0 = h(V ,W ) s s s s or e q u i v a l e n t l y , U = a,V / [ a 2 + a,V ), (4.1-1) s s s W = a 2V / [ a , + a 2V ]. (4.1-2) s s s Moreover, (U ,V ,W ) s a t i s f i e s the c o n s e r v a t i o n equation (4.5), s s s U / ( a , 6 ) + V + (a,/c)W = 5, (4.1-3) s s s where $ = U 0 /(a,5) + V 0 + (a,/e)W 0 + j3(t*)/L (4.1-4) and /3(t*) i s d e f i n e d by equation (4.4). Using equations (4.1 - 1)-(4.1-3) we can prove the f o l l o w i n g p r o p o s i t i o n : P r o p o s i t i o n II : For a l l p o s i t i v e parameter values a 1 f a 2 , 5 , and e, the s t e a d y - s t a t e s o l u t i o n to system (P) i s unique. 26 Proof : Observe that equations (4.1 -1)-(4.1-2) d e f i n e a m o n o t o n i c a l l y i n c r e a s i n g curve i n a 3-dimensional space with V as a parameter s ranging from 0 to i n f i n i t y (see F i g u r e 2 ) . Thus the curve i n t e r s e c t s the plane d e f i n e d by equation (4.1-3) a t p r e c i s e l y one p o i n t , (U , V ,W ) , and t h i s completes the p r o o f , s s s V F i g u r e 2 : Geometrical i l l u s t r a t i o n of the uniqueness of the s t e a d y - s t a t e s o l u t i o n . In order t o f i n d the s t e a d y - s t a t e v a l u e (U , V ,W ) , we s s s s u b s t i t u t e equations (4.1 -1)-(4.1-2) i n t o (4.1-3) and we see t h a t V must s a t i s f y the c u b i c equation s V 3 + / c , V 2 + K 2 V - S = 0, (4.1-5) s s s 27 where «, = a,/e + 1/(83,) + ( a , / a 2 + a 2/a,) - $, K2 = 1 + 1/(8a 2) + a 2/e - ^(a,/a2 + a 2 / a , ) . By P r o p o s i t i o n I I , (4.1-5) has a unique p o s i t i v e r e a l r o o t . Using formulas from a standard t a b l e [25], we o b t a i n V = fc,/3 + A, + A 2 , (4.1-6) s where A, 3 = -p/2 + v/p2/4 + q 3/27, A 2 3 = -p/2 - /p 2/4 + q 3/27, and where 1 p = [ 2 K , 3 - 9K,K2 + 27$]/27, q = [ 3 K 2 - K , 2 ] / 3 . The v a l u e s of U and W can now be obtained from equations s s (4.1 - 1)-(4.1-2) . For the set of parameter values shown in Table 2, the s t e a d y - s t a t e value i s U = 0.00711, V = 0.04525, W = 0.22249. s s s In p a s s i n g we remark that due to the complexity of the c o e f f i c i e n t s K , and K2 , i t i s d i f f i c u l t to prove P r o p o s i t i o n II d i r e c t l y from using equation (4.1-5) and t h i s i l l u s t r a t e s the importance of g e o m e t r i c a l ideas i n mathematics. We s h a l l now turn to study the s t a b i l i t y of the s t e a d y - s t a t e s o l u t i o n . 1Note that s i n c e (4.1-5) has a unique p o s i t i v e r e a l r o o t , the e x p r e s s i o n p 2/4 + q 3/27 > 0. 28 4.2 ASYMPTOTIC STABILITY OF THE STEADY-STATE We s h a l l begin by g i v i n g a formal d e f i n i t i o n of s t a b i l i t y . Let y (x) be a s t e a d y - s t a t e s o l u t i o n to a system of PDEs of s _ type (1.2) with i n i t i a l c o n d i t i o n y 0 ( x , t 0 ) a n <3 a p p r o p r i a t e boundary c o n d i t i o n s . Let |«| denote the usual E u c l i d e a n norm. The s t e a d y - s t a t e y (x) i s s a i d t o be s t a b l e i f given any e > 0, s there e x i s t s a 6 such that | y 0 ( x , t 0 ) - y ( x ) | < 6 => |y(x,t) - y ( x ) | < e, s s fo r a l l t > t 0 . I f i n a d d i t i o n |y(x,t) - y (x)| -> 0 as t -> °°, s then the s t e a d y - s t a t e y (x) i s a s y m p t o t i c a l l y s t a b l e . s I n t r i n s i c a l l y , i t i s e a s i e r to prove i n s t a b i l i t y than i t i s to e s t a b l i s h s t a b i l i t y . For to e s t a b l i s h s t a b i l i t y we have to show that every s u f f i c i e n t l y s m a l l p e r t u r b a t i o n d i e s away, whereas i n s t a b i l i t y w i l l be demonstrated i f we can show that one s u f f i c i e n t l y small p e r t u r b a t i o n away from the s t e a d y - s t a t e w i l l grow with time. S u f f i c i e n t c o n d i t i o n s f o r i n s t a b i l i t y are given by Amundson [3] and Jackson [17]. In t h i s s e c t i o n , however, we s h a l l use the l i n e a r i z a t i o n t e c h n i q u e 1 to e s t a b l i s h the asymptotic s t a b i l i t y of (U ,V ,W ) s s s ' J u s t i f i c a t i o n of the l i n e a r i z a t i o n technique f o r e s t a b l i s h i n g s t a b i l i t y i s d i s c u s s e d by K r a s n o s e l ' s k i i [19]. 29 found i n the p r e v i o u s s e c t i o n . The l i n e a r i z e d form of the governing equations i s 9u/9t = D 1 3 9 2 u / 9 x 2 - au + a , bv, 9v/9t = 9 2v/9x 2 + au/(a , 6 ) - (b/5 + c/e)v + (da,/e)w f 9w/9t = D 2 39 2w/9x 2 + (c/a,)v - dw, (4.2-1) where a = a,V + a 2 , b = ( l - U ) , s s c = (1 - W ), d = (V /a, + 1/a 2), s s and where (u,v,w) are d e v i a t i o n s from (U ,V ,W ), i . e . , s s s u = U - U , v = V - V , w = W - W . s s s Moreover f o r every given s t e a d y - s t a t e value s a t i s f y i n g (4.5), i t i s easy to deduce from (4.3) that the d e v i a t i o n must s a t i s f y the i n t e g r a l equation L / [ u ( x , t ) / ( a , 5 ) + v ( x , t ) + (a,/e)w(x,t)] dx = 0. (4.2-2) 0 Since the system of equations (4.2-1) i s l i n e a r and has constant c o e f f i c i e n t s , i t s s o l u t i o n can be found by s e p a r a t i o n of v a r i a b l e s . Moreover, f o r t > t * , the n o n l i n e a r boundary c o n d i t i o n s (3.3) reduce to the homogeneous Neumann type and hence we w r i t e u(x,t) = Z u (t)cos(w x ) , v ( x , t ) = I v ( t ) c o s ( u x ) , n=0 n n n=0 n n w(x,t) = Z w (t)cos(w x ) , (4.2-3) n = 0 n n 30 where u> = n 7 r / L , n = 0 , 1 , 2 , . . . n S u b s t i t u t i n g equation (4.2-3) i n t o (4.2-1) we o b t a i n a system of o r d i n a r y - d i f f e r e n t i a l equations z ' (t) = A z ( t ) , (4.2-4) n n n where z (t) = (u ( t ) , v (t),w ( t ) ) T and n n n n -(D, 3o 2 + a) a,b 0 n A = a/Ua,) -(CJ2 + b/6 + c/e) da,/e n n 0 c/a, - (D 2 3CJ2 + d) We s h a l l f i r s t e s t a b l i s h the f o l l o w i n g property of matrix A , n Property 4 : For a l l p o s i t i v e parameter values a,, a 2 , 6, and e, a l l the eigenvalues of matrix A have negative r e a l p a r t s f o r n > 0 . When n n = 0 , one of the eigenvalues i s zero and the other two eig e n v a l u e s have negative r e a l p a r t s . Proof : To f i n d the eigenvalues of A we set n f(X) = | XI - A | = X 3 + aX 2 + /3X + 7 = 0 , ( 4 . 2 - 5 ) n where I i s the 3x3 i d e n t i t y matrix and 31 a = (1 + D 1 3 + D 2 3)CJ 2 + (a + b/5 + c / e + d) , n 0 = ( D 1 3 + D 1 3 D 2 3 + D 2 3 ) u 4 n + [(b/5 + c / e + d ) D 1 3 + (a + b/6 + c / e ) D 2 3 + (a + d) ]w2 n + [ac/e + ad + bd/6], 7 = D 1 3D 2 3w 6 + [ d D 1 3 + (b/6 + c / e ) D 1 3 D 2 3 + a D 2 3 K n n + [ ( b d / 6 ) D 1 3 + ad + (ac / e )D 2 3 ]co2 . n Observe that s i n c e a, b, c, and d are a l l p o s i t i v e , so are a, j3, and 7. When n = 0, we have 7 = 0 and the r e q u i r e d r e s u l t f o l l o w s immediately s i n c e a and j3 are p o s i t i v e . For n > 0, observe that f ( 0 ) = 7 > 0, and f ( - a ) = 7 - a/3 < 0. Thus there e x i s t s a negative r e a l number X,, -a < X, < 0, such that f ( X , ) = 0. D i v i d i n g equation ( 4 . 2 - 5 ) by i t s f a c t o r ( X - X , ) , we see that the other two roots must s a t i s f y the q u a d r a t i c equation X 2 + ( a + X , ) X + [ 0 + X T U + X , ) ] = 0. Now s i n c e a + X, > 0, the r o o t s to the q u a d r a t i c must have negative r e a l p a r t s and t h i s completes the p r o o f . According to the standard theory of l i n e a r s t a b i l i t y [ 7], Property 4 s t a t e s that f o r n > 0, the system of o r d i n a r y d i f f e r e n t i a l equations ( 4 . 2 - 4 ) i s s t a b l e and the s o l u t i o n z (t) = exp(A t) -> 0 as t -> °°. For n = 0, ( 4 . 2 - 4 ) i s n e u t r a l l y n n 32 s t a b l e and the s o l u t i o n z 0 ( t ) = e x p ( A 0 t ) i s bounded. Using the above r e s u l t we can prove the f o l l o w i n g p r o p o s i t i o n about the s t a b i l i t y of system (P). P r o p o s i t i o n III : For a l l p o s i t i v e parameter val u e s a,, a 2 , 5, and e, the unique s t e a d y - s t a t e s o l u t i o n (U ,V ,W ) i s a s y m p t o t i c a l l y s t a b l e . s s s Proof : It s u f f i c e s to show that the d e v i a t i o n (u,v,w) tends to zero in the s t e a d y - s t a t e . By Property 4 the s t e a d y - s t a t e value can be obtained by merely s o l v i n g z 0 ' ( t ) = A 0 z 0 ( t ) and l e t t i n g t —> °°. The s o l u t i o n i s u(x,<=°) = u 0 ( » ) = ( a , b / a ) M , v(x,°°) = v 0 (°°) = n, w(x,o°) = w 0 ( » ) = ( c/da,)n, (4.2-6) where M i s some constant depending on the i n i t i a l c o n d i t i o n as w e l l as the parameters a,, a 2 , 6, and e. Since u(x,°°), v(x,°°) and w(x,«) are independent of x, the i n t e g r a l equation (4.2-2) reduces to the a l g e b r a i c equation u 0(»)/(a,5) + v0(=°) + (a,/e)w 0(») = 0. (4.2-7) Combining equations (4.2-6) and (4.2-7) we deduce that 33 u(x,°°) = v(x,°°) = w(x,°°) = 0 and t h i s completes the pro o f . 34 5. SOME FEATURES OF THE TRANSIENT BEHAVIOUR One of our main o b j e c t i v e s i n f o r m u l a t i n g the mathematical model (P) i s to provide q u a l i t a t i v e and q u a n t i t a t i v e i n f o r m a t i o n about the behaviour of the i n i t i a l t r a n s i e n t changes of C a 2 + . In t h i s chapter we s h a l l approximate the t r a n s i e n t s o l u t i o n v i a a reg u l a r p e r t u r b a t i o n approach. Moreover, we s h a l l use a comparison theorem to o b t a i n a n a l y t i c a l upper and lower bounds on the exact s o l u t i o n . We mentioned i n S e c t i o n 4.1 that d u r i n g the i n i t i a l t r a n s i e n t s t a t e , the d i f f u s i o n r a t e of system (P) i s very r a p i d compared with i t s r e a c t i o n r a t e . Thus before c a r r y i n g out the p e r t u r b a t i o n a n a l y s i s , we s h a l l r e s c a l e our v a r i a b l e s so that the d i f f u s i o n terms are of dominant order. To t h i s end we introduce the s t r e t c h i n g v a r i a b l e s 1 H = ,/(5/e 2)x, r = ( 6 / e 2 ) t , UU,T) = U / ( e / 5 ) 2 , V U , T ) = V/U/5) , w(£,r) = W/U/5). In (£,T) space, system (P) becomes 9u/9r = D 1 3 9 2 u / 9 £ 2 + (e/6)[5a,v] - ( e / 5 ) 2 [ 5 a 2 u ] + (e/8)3[6a,uv], 9V/9T = 9 2v/9£ 2 + (e/5)[a,w/a 2 - v] + ( e / 5 ) 2 [ a 2 u / a , - V(1-W)] + U / 5 ) 3 [ u v ] , 1The new dependent v a r i a b l e s are chosen such that the i n i t i a l c o n d i t i o n s are of 0 ( 1 ) . 35 9w/9r = D 2 3 9 2 w / 9 £ 2 + ( e / 5 ) 2 [ 6 ( v / a , - w/a 2)] - (e/5) 3[6wv/a,], (5.1) with i n i t i a l c o n d i t i o n s u(£,0) = u 0 , v(£,0) = v 0 , w(£,0) = w0, (5.2) where u 0 = U 0 / ( e / 5 ) 2 , v 0 = V 0 / ( e / 5 ) , w0 = W 0 / ( e/5). The boundary c o n d i t i o n s are 9U(0,T)/9$ = 9u(L',r)/9£ = 9w(0,r)/9£ = 9w(L',r)/9$ = 0, 9v(L*,r)/9£ = 0, 9v(0,r)/9£ = A' + B ' l n v ( 0 , T ) , (5.3) where A* = /6(A + Bin ( e / 5 ) ) , B' = /SB, L' = i / ( 5 / e z ) L . I t i s c l e a r from (5.1) that (e/6) i s a n a t u r a l choice as a p e r t u r b a t i o n parameter. For the set of parameter values shown i n Table 2, e/5 has a value of 0.0375. In order to o b t a i n an approximate t r a n s i e n t s o l u t i o n to system (P), we c o n s i d e r a r e g u l a r expansion of the form UU,T) = U°U,T) + ( e / 6 ) u 1 U , T ) + ( e/6) 2u 2(£,r) + ... v U , r ) = V ° ( { , T ) + ( e / 6 ) v 1 U , r ) + ( e/5 ) 2 v 2 ( % , r ) + ... W($,T) = W°U,T) + (e/6)w'({,r) + ( e/6) 2w 2 ( £ , r) + ... (5.4) The n o n l i n e a r boundary c o n d i t i o n on v(£,r) becomes A' + B'lnv(0,T) = A' + B'lnv°(0,r) + (e/5)[B'/v°(0,r)]v 1(o,r) + 0 ( ( e / 6 ) 2 , ) . 36 Our o b j e c t i v e i s to o b t a i n the f i r s t - o r d e r c o r r e c t i o n terms in the expansion (5.4). S u b s t i t u t i n g equation (5.4) i n t o (5.1) and c o l l e c t i n g terms of the same order g i v e s 0((e/6)°) : 3u°/3r = D 1 3 3 2 u ° / 3 £ 2 , 3v°/9r = 3 2 v ° / a $ 2 , 3w°/3r = D 2 3 3 2 w ° / 9 £ 2 , (5.5) with i n i t i a l and boundary c o n d i t i o n s u ° U , 0 ) = u 0 , v ° U , 0 ) = v 0 , w ° U,0) = w0, 3U°(0,T)/3$ = 3u°(L',r)/3£ = 3w°(0,r)/3£ = 3w°(L' ,r)/3 £ = 0, 3V°(L',T)/3£ = 0, 3V°(0,T)/3£ = A' + B'lnv°(0,r). (5.6) Thus the zeroth-order e q u a t i o n s 1 to our n o n l i n e a r system are l i n e a r uncoupled heat equations with n o n l i n e a r boundary c o n d i t i o n on v 0 ( £ , r ) . In view of the i n i t i a l and boundary c o n d i t i o n s , the s o l u t i o n s to U°(£,T) and w°(£,r) can be obtained by i n s p e c t i o n as u ° U , r ) = u 0 , w ° U , r ) = w0. (5.7) However, because of the n o n l i n e a r boundary c o n d i t i o n , the exact s o l u t i o n to v°(£,r) cannot be obtained a n a l y t i c a l l y . Moreover, l i n e a r i z a t i o n of the boundary c o n d i t i o n i s not f r u i t f u l s i n c e v°(£,r) changes r a p i d l y near T = 0. Thus, v°(£,r) i s c a l c u l a t e d n u m e r i c a l l y using an a l g o r i t h m c a l l e d PDECOL (see 10n p h y s i c a l grounds, the ze r o t h - o r d e r s o l u t i o n corresponds to an approximate s o l u t i o n of the system when the d i f f u s i o n r a t e i s r a p i d r e l a t i v e to the r e a c t i o n r a t e . 37 Appendix A f o r a d e t a i l d i s c u s s i o n ) . We s h a l l proceed to ob t a i n c o r r e c t i o n terms f o r u, v, and w. The f i r s t - o r d e r governing equations are 0 ( ( e / S ) ) : 3u 1/3r = D 13 3 2 u V 3 £ 2 + 6a,v°, (5.8a) 3V 1/3T = 3 2 v 1 / 3 £ 2 + (a,w°/a 2 - v ° ) , (5.8b) 3w'/3r = D 2 3 3 2 w 1 / 3 £ 2 , (5.8c) with i n i t i a l and boundary c o n d i t i o n s u M ^ O ) = 0, v ' ( L O ) = 0, w ' ( L O ) = 0, 3U 1(0,T)/3£ = 3U 1(L*,T)/3£ = 9 W ' ( 0 , T ) / 3 J = B w M L ' . r J / B J = 0, 3 v 1 ( L \ r ) / 3 £ = 0, 3v 1 ( 0 , r ) / 3 £ = [ B ' /v° ( 0 , r ) ] v 1 ( 0 , r ) . (5.9) An i n s p e c t i o n of the i n i t i a l and boundary c o n d i t i o n s on w1 r e v e a l s that w1 U , r ) = 0. (5.10) Thus i n order to obt a i n a c o r r e c t i o n term f o r W(£,T), we must c o n s i d e r the 0 ( ( e/6) 2) equation and we have 0 ( ( e/6) 2) : 3W2/3T = D 2 33 2w 2/3$ 2 + 6 [ v % , - w ° / a 2 ] , (5.11) with i n i t i a l and boundary c o n d i t i o n s w 2(£,0) = 0, 3w 2(0,r)/3£ = 0 = 3w 2(L',T)/3£. (5.12) Equations (5.8a) and (5.11) are p r e c i s e l y of the form given i n Appendix B with a = 0. The s o l u t i o n s are u 1 U , r ) = £ u (T)COSX ( £ - L ' ) , (5.13a) n = 0 n n 38 00 w 2 U , r ) = I w ( T ) C O S X U-L') - ( 8 w 0 / a 2 ) r , (5.13b) n = 0 n n where X = n.7r/L' , n r u ( T ) = [5a,] / q ( s ) e x p ( - D , 3 X 2 ( r - s ) ) d s , n O n n r w (T) = [ 6 / a J J q ( s ) e x p ( - D 2 3 X 2 ( r - s ) ) d s , n O n n fo r n = 0,1,2,... and where L' q 0 ( s ) = (1/ L ' ) / v ° U , s ) d £ , 0 L' q (s) = ( 2 / L ' ) ; v°(£,s)cosX (£-L')d£, n > 1. n 0 n Since v°(£,r) i s known from the zeroth-order c a l c u l a t i o n , q (s) n can be obtained by d i r e c t numerical i n t e g r a t i o n . We s h a l l now turn to the d i s c u s s i o n of v 1 ( £ , r ) . By Lemma 1 of Appendix B, provided that |A[B'/v°(0,r)]/[B'/(v°(0,r))]| « 1, or e q u i v a l e n t l y |[v°(0,r+Ar) - V°(0,T)]/v°(0,r)| « 1, (*) then v 1(£,r) i s approximately given by 00 v 1 U , r ) = I v (r)cosX U - L ' ) , (5.14) n=1 n n where X s a t i s f i e s the t r a n s c e n d e n t a l equation n X sinX L' - [B'/v°(0,r)]cosX L' = 0, (5.15) n n n 39 and wi th and where T v (T) = / p ( s ) e x p ( - X 2 ( r - s ) ) d s , n O n n L 1 p (s) = [J ( a i w 0 / a 2 - v ° U , s ) ) c o s X U-L')d£]/c , n 0 n n c = { L 1 + [ B ' / V ° ( 0 , T ) ] c o s 2 X L'/X 2}/2, n n n for n = 1,2,. In order to ensure that c o n d i t i o n (*) holds we d i v i d e the i n t e r v a l [0,T] i n t o m s u b i n t e r v a l s 0 = T < T < < T = T , 0 1 m so that T m k v ( T) = Z J p ( s ) e x p ( - X 2 ( r - s ) ) d s . n k=1 r n n k-1 For each s u b i n t e r v a l r < r < T , we approximate v°(0,r) by k-1 k V°(0,T ) and generate a new set of {X } v i a equation (5.15). k-1 n Thus by choosing m s u f f i c i e n t l y l a r g e , c o n d i t i o n ( * ) i s s a t i s f i e d . In summary the approximate t r a n s i e n t s o l u t i o n s obtained from the p e r t u r b a t i o n method are —" u U , r ) = u 0 + ( e / 5 ) u 1 U , r ) + 0 ( ( e / 6 ) 2 ) , v U , r ) = v ° U , r ) + ( e/5 ) v ' ( { , r ) + 0 ( ( e / 5 ) 2 ) , w(^,r) = w0 + ( e/5) 2w 2(^,r) + 0 ( ( e / 5 ) 3 ) . (5.16) 40 In terms of U, V and W, the p e r t u r b a t i o n s o l u t i o n s are U ( $ , T ) = Uo + ( e / 6 ) 3 u 1 U , r ) + 0 ( ( e / 5 ) M , (5.17a) V ( £ , r ) = V ° U , T ) + (e/S) 2v'(£,r) + 0 ( U / 5 ) 3 ) , (5.17b) W(£,r) = W0 + (e/8) 3w 2(£,r) + 0 ( ( e / 6 ) a ) , (5.17c) where V°(£,r) = (e/5)v°(£,7). Thus our p e r t u r b a t i o n s o l u t i o n s are good up to 0 ( ( e / 6 ) 2 ) fo r V and 0 ( ( e / 6 ) 3 ) f o r U and W. For v a l u e s 1 of t ^ 10~ 3 (t = r ( e 2 / 8 ) ) , equations (5.17a-c) provide good approximations to the i n i t i a l t r a n s i e n t behaviour of system (P). In F i g u r e s (3a-b) we compare the exact s o l u t i o n of V ( £ , T ) with the p e r t u r b a t i o n s o l u t i o n s f o r t = 10" 3 and 10" 2. We have omitted p l o t t i n g the U(£,r) and W(£,r) s o l u t i o n s s i n c e f o r values of t < 10" 2, the d e v i a t i o n s from the i n i t i a l s t a t e , U 0 and W0, are n e g l i g i b l e . Our p e r t u r b a t i o n a n a l y s i s i s u n s a t i s f a c t o r y in the sense that we are not able to o b t a i n an a n a l y t i c a l s o l u t i o n 2 f o r V ° ( £ , r ) . Moreover for values of t > 10~ 3, the p e r t u r b a t i o n s o l u t i o n s (5.17a-c) begin to d e v i a t e from the exact s o l u t i o n s and they no longer provide any q u a l i t a t i v e i n f o r m a t i o n of the true s o l u t i o n s . Thus one would at l e a s t l i k e to f i n d good a n a l y t i c a l bounds which w i l l g ive us q u a n t i t a t i v e i n f o r m a t i o n on the exact s o l u t i o n s . The f o l l o w i n g v e r s i o n of the comparison theorem [10] f o r p a r a b o l i c equations s u i t s our s i t u a t i o n . 1The value t = 10" 3 corresponds to 1 Msec in p h y s i c a l u n i t s . 2However the computing time f o r s o l v i n g V°(£,r) i s l e s s than h a l f the computing time f o r s o l v i n g the f u l l system of PDEs. F i g u r e 3a-b : Comparison between the exact s o l u t i o n of V with the p e r t u r b a t i o n s o l u t i o n f o r t = 0.001, 0.01. 42 Comparison Theorem Consider the system of p a r a b o l i c equations of type (1.2) 3v/3t = D3 2v/3x 2 + g'(v) , 3u/3t = D3 2u/3x 2 + F ( u ) , d e f i n e d f o r x e I = [a,b], t > 0 with p r e s c r i b e d i n i t i a l c o n d i t i o n s and Neumann boundary c o n d i t i o n s . Assume i ) F(u) i s quasi-monotone nondecreasing in u, i i ) G(-) > F ( - ) , i i i ) v(x,0) > u(x,0), x e I = [a, b ] , iv ) 3v(x,t)/3n > 3u(x,t)/3n, x e 31, t > 0, where 3/3n denotes outward d i r e c t i o n a l d e r i v a t i v e on 31. Then v ( x , t ) > u ( x , t ) , f o r x e I, t > 0, for any continuous s o l u t i o n s u(x,t) and v ( x , t ) . Thus p r o v i d e d that c o n d i t i o n ( i ) i s s a t i s f i e d , the above theorem a l l o w s us to c a l c u l a t e an upper bound f o r u by r e p l a c i n g the k i n e t i c term, the i n i t i a l c o n d i t i o n , and the normal d e r i v a t i v e s by any convenient d i f f e r e n t i a b l e f u n c t i o n s s u b j e c t to c o n d i t i o n s ( i i ) - ( i v ) . We s h a l l use the above theorem to f i n d a n a l y t i c a l bounds f o r U, V, and W. In p a r t i c u l a r we are i n t e r e s t e d i n f i n d i n g good bounds d u r i n g the i n i t i a l t r a n s i e n t s t a t e . By Property 3 of Chapter 4, we know that U 0, V 0, and W0 are lower bounds f o r U, V, and W, r e s p e c t i v e l y . Moreover our p e r t u r b a t i o n s o l u t i o n s i n d i c a t e 43 that f o r small t , they are good lower bounds. Thus we s h a l l merely c o n c e n t r a t e on f i n d i n g a n a l y t i c a l upper bounds. The f o l l o w i n g upper bounds f o r f, g, and h are easy to e s t a b l i s h : f(U,V) < a,(1-U 0)(V -V 0) = f, (5.19a) m _ g(U,V,W) < a 2/(5a,) + a , / ( e a 2 ) = g, (5.19c) h(V,W) < (1-W 0)(V - V 0 ) / a , = h, (5.19b) m where V i s a f i n i t e upper bound on V to be determined, m We s h a l l t e m p o r a r i l y assume that V i s known and prove the m f o l l o w i n g p r o p o s i t i o n : P r o p o s i t i o n V : A pointwise upper bound (U,V,W) f o r (U,V,W) i s given by U(x,t) = U 0 + f t , (5.20a) W(x,t) = W0 + ht, (5.20b) V(x,t) = V 0 + gt + V * ( x , t ) , (5.20c) where OO OO V * ( x , t ) = 2/3/t { Z ( l/ v/7i :)exp(-(x-2kL) 2 / 4 t ) - Z P[ ( 2kL-x ) /2 f t ] k=-« k=1 OO - Z P[(x+2kL)/2/t]}, (5.21) k = 0 where /3 = -(A + B l n V 0 ) , P(z) = z e r f c ( z ) and e r f c ( x ) i s the complementary e r r o r f u n c t i o n , OO e r f c ( x ) = 2/V7T / e x p ( - a 2 ) d a . x 44 Proof : By Property 1 of Chapter 4, ¥ - (f,g,h) i s quasi-monotone nondecreasing i n u = (U,V,W). Thus c o n d i t i o n ( i ) of our comparison theorem i s s a t i s f i e d and an upper bound to system (P) can be ob t a i n e d by s o l v i n g the system of equations 3U/3t = D 1 33 2U/3x 2 + I , (5.22a) 3V/3t = 3 2V/3x 2 + g, (5.22b) 3W/3t = D 2 33 2W/3x 2 + h, (5.22c) with i n i t i a l c o n d i t i o n s on 0 < x < L, U(x,0) = U 0, V(x,0) = V 0, W(x,0) = W0, (5.23) and boundary c o n d i t i o n s f o r t > 0, 3U/3x(0,t) = 0 = 3U/3x(L,t), 3W/3x(0,t) = 0 = 3W/3x(L,t), 3V/3x(0,t) = -/3, 3V/3x(L,t) = 0. (5.24) Observe that we have kept the same i n i t i a l c o n d i t i o n s and have simply r e p l a c e d the k i n e t i c terms by t h e i r r e s p e c t i v e upper bounds. Now s i n c e 3V/3x(0,t) = -j3 = A + Bln V 0 < A + BlnV(0,t) = 3V/3x(0,t), a l l c o n d i t i o n s of the theorem are s a t i s f i e d . The s o l u t i o n s f o r U and W can be obtained by i n s p e c t i o n as U(x,t) = U 0 + f t , W(x,t) = W0 + ht. To so l v e f o r V(x,t) we c o n s i d e r the t r a n s f o r m a t i o n V * ( x , t ) = V - V 0 - gt. (5.25) 45 S u b s t i t u t i n g (5.25) i n t o (5.21b) we have 3V*/3t = 3 2V*/3x 2, (5.26) subj e c t to i n i t i a l c o n d i t i o n on 0 < x < L, . V*(x,0) = 0, (5.27) and boundary c o n d i t i o n s f o r t > 0, 9V*/3x(0,t) = -/3, 3V*/3x(L,t) = 0. (5.28) We s h a l l s o l v e (5.26) using the Laplace transform. Taking the transform on both s i d e s of (5.26) and using i n i t i a l c o n d i t i o n (5.27) we have d 2v*/dx 2 - sv* = 0, (5.29) where v*(x,s) = £{V*(x,t)}. The boundary c o n d i t i o n s are dv*(0,s)/dx = -p/s, ' dv*(L,s)/dx = 0. (5.30) S o l v i n g (5.29) sub j e c t to c o n d i t i o n s (5.30) g i v e s v*(x,s) = (p/^/J7) [coshv/s(x-L)/sinhv/sL]. (5.31) Since we are i n t e r e s t e d i n the s o l u t i o n f o r small t , we expand (5.31) f o r l a r g e s and we get CO OO v*(x,s) = (p/\/sT){ L exp ( - V l(2kL-x) ) + Z exp(-/s (x + 2kL))}. k=1 k=0 (5.32) Using a standard transform t a b l e [25], we can i n v e r t (5.32) term by term and we o b t a i n (5.21) as d e s i r e d . 46 An i n t e r e s t i n g f e a t u r e of the i n f i n i t e s e r i e s r e p r e s e n t a t i o n of V * ( x , t ) i s t h a t i t converges r a p i d l y f o r s m a l l t . To see t h i s r e c a l l t h a t the a s y m p t o t i c e x p r e s s i o n f o r e r f c ( x ) i s e r f c ( x ) = ( l / v / 7 r ) e x p ( - x 2 ) [ 1/x - 1/2X 3 + 0 ( 1 / X 5 ) ] as x -> +». Us i n g the above e x p a n s i o n i t i s easy t o deduce t h a t CO V * ( x , t ) = 0/17* ( Z e x p ( - ( x - 2 k L ) 2 / 4 t ) / [ ( x - 2 k L ) 2 / 4 t ] k = - c o + 0 ( [ ( 2 k L - x ) / 2 / t ] " 5 ) } as [ ( 2 k L - x ) / 2 / t ] -> +». (5.33) For v a l u e s of t < 0.1, i t i s o n l y n e c e s s a r y t o c a l c u l a t e the k = 0,±1 terms t o o b t a i n s u f f i c i e n t a c c u r a c y f o r V * ( x , t ) . I t i s c l e a r from (5.33) t h a t f o r s m a l l t , t h e c o n t r i b u t i o n of V* i n (5.20c) i s s m a l l 1 compared w i t h V 0 + g t . Thus f o r v a l u e s of t < 1, we s h a l l n e g l e c t the V* term i n (5.20c) and we have V ( t ) = V ( x , t ) = V 0 + g t , x e [ 0 , L ] , 0 < t < 1. m Thus f o r s m a l l t , the bounds (U,V,W) are independent of x. U s i n g the s e t of parameter v a l u e s g i v e n i n Table 1, we summarize i n T a b l e 3 the t r u e maximum v a l u e s of (U,V,W) which occur a t x = 0 t o g e t h e r w i t h t h e i r r e s p e c t i v e upper bounds f o r t = 0.001, 0.01, and 0.1. —~ 'Numerical c a l c u l a t i o n s i n d i c a t e t h a t f o r v a l u e s of t < 1.0, the c o n t r i b u t i o n from V* i s n e g l i g i b l e . 47 t = 0.001 t = 0.01 t = 0. 1 Max U(x,t) x U 0.00250 0.00250 0.00258 0.00374 0.00314 0. 12719 Max V(x,t) X V 0.07776 0.09488 0.13423 0.80651 0. 14558 7.92283 Max W(x,t) X W 0.09111 0.09136 0.09514 0.13638 0. 12851 4.63752 Table 3 : Summary of the true maximum values of (U,V,W) together with t h e i r r e s p e c t i v e upper bounds. It i s c l e a r from Table 3 that as t i n c r e a s e s the bounds d e v i a t e f u r t h e r and f u r t h e r away from the exact s o l u t i o n s and hence they are not very u s e f u l f o r l a r g e values of t . T h i s i s due mainly to the rough estimate we used i n c a l c u l a t i n g f, g, and h in equations ( 5 . l 9 a - c ) . Our present a n a l y s i s , however, p r o v i d e s us with an idea of the r e l a t i v e maximum amplitude of U, V, and W. We s h a l l now turn to d i s c u s s the numerical r e s u l t s of system (P). 48 6. NUMERICAL SOLUTION In the area of c h e m i c a l l y r e a c t i n g systems, n o n l i n e a r i t y i s more the r u l e than the e x c e p t i o n , and a c c o r d i n g l y a n a l y t i c a l s o l u t i o n s are r a r e . We were " f o r t u n a t e " to o b t a i n a n a l y t i c a l r e s u l t s f o r the s t e a d y - s t a t e and the i n i t i a l t r a n s i e n t behaviour of system ( P ) . However f o r int e r m e d i a t e values of t , a n a l y t i c a l r e s u l t s f o r system (P) are beyond our reach and we must r e s o r t to numerical techniques. In order to give a complete d e s c r i p t i o n of the c a l c i u m dynamics, we use the parameter val u e s i n Table 1 and solve the f u l l system of equations (3.2)-(3.4) using PDECOL. F i g u r e s (4a-c) show the c o n c e n t r a t i o n p r o f i l e s of U, V, and W .([CaArz], [ C a 2 + ] , and [CaB]) i n t h e i r o r i g i n a l p h y s i c a l u n i t s d u r i n g a 0.1 second (t * ) l o a d i n g p e r i o d . I t i s c l e a r from the f i g u r e s that the c o n c e n t r a t i o n s of CaArz, C a 2 + , and CaB r i s e to t h e i r r e s p e c t i v e maxima 1 at t = t * and as C a 2 + i n f l u x c o n t i n u e s , the f r o n t of CaB and CaArz moves away from the s u r f a c e membrane (or i n t o the c e l l ) . In F i g u r e s (5a~c) we show that a f t e r the pulse stops (t > t * ) , d i f f u s i o n proceeds c a r r y i n g C a 2 + and the complexes f u r t h e r i n t o the c e l l which permits a s p a t i a l r e d i s t r i b u t i o n of C a 2 + i n s i d e the cytoplasm. The c o n c e n t r a t i o n s 1At the same time, the i n t r i n s i c b u f f e r c o n c e n t r a t i o n s Arz and B decrease a c c o r d i n g to equations (2.2a) and (2.2b), r e s p e c t i v e l y . 49 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 x (ym) F i g u r e 4a : C o n c e n t r a t i o n p r o f i l e s of CaArz d u r i n g a 0.1 second l o a d i n g p e r i o d . 50 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 X ( l im) F i g u r e 4b : C o n c e n t r a t i o n p r o f i l e s of C a 2 + d u r i n g a 0.1 second l o a d i n g p e r i o d . 51 F i g u r e 4c : C o n c e n t r a t i o n p r o f i l e s of CaB d u r i n g a 0.1 second l o a d i n g p e r i o d . 52 F i g u r e 5a : C o n c e n t r a t i o n p r o f i l e s of CaArz a f t e r i n f l u x . 53 [Ca] vM/1 x (ym) F i g u r e 5b : C o n c e n t r a t i o n p r o f i l e s of C a 2 + a f t e r i n f l u x . 54 [CaB] yM/1 \t = 100 msec 14.0 x (ym) F i g u r e 5c : C o n c e n t r a t i o n p r o f i l e s of CaB a f t e r i n f l u x . 55 of the C a 2 + and the complexes s t e a d i l y approach t h e i r r e s p e c t i v e e q u i l i b r i u m v a l u e s and the recovery time i s approximately 8 seconds. Using the set of parameter values given i n Table 1, we summarize i n Table 4 the i n i t i a l , maximum, and s t e a d y - s t a t e c o n c e n t r a t i o n l e v e l s of CaArz, C a 2 + , and CaB. t sees [CaArz] MM/1 [ C a 2 + ] MM/1 [CaB] MM/1 I n i t i a l s t a t e 0.0 0.7482 0.1000 18.1818 Maximum value 0.1 49.0101 10.0776 181.1472 Steady-state value 8.0 2.1413 0.2876 44.6679 Table 4 : Summary of the i n i t i a l , maximum, and s t e a d y - s t a t e c o n c e n t r a t i o n l e v e l s of CaArz, C a 2 + , and CaB. F i n a l l y we remark that our numerical r e s u l t s are i d e n t i c a l to those given by Connor and Nikolakopoulou. (See F i g u r e s 5A-B and F i g u r e 7 of t h e i r paper [13].) As mentioned i n t h e i r paper, the numerical r e s u l t s do not only give a good q u a l i t a t i v e account of the c a l c i u m dynamics, but q u a n t i t a t i v e l y , they are w i t h i n 50% of the experimental v a l u e s . P o s s i b l e changes i n the model parameters which give b e t t e r experimental-computational agreement are a l s o d i s c u s s e d i n t h e i r paper. 56 7. CONCLUSION In t h i s t h e s i s we analyse a d i f f u s i o n model which a r i s e s in neurobiology d e s c r i b i n g the dynamics of i n t r a c e l l u l a r c a l c i u m ion c o n c e n t r a t i o n changes due to d i f f u s i o n and to b u f f e r i n g i n nerve cytoplasm. The model was o r i g i n a l l y formulated by Connor and Nikolakopoulou [13] i n an e f f o r t to answer some q u a n t i t a t i v e q u e s t i o n s i n c l u d i n g the s p a t i a l d i s t r i b u t i o n of ca l c i u m w i t h i n the cytoplasm and the i n c r e a s e i n i n t r a c e l l u l a r [ C a 2 + ] that a given i n f l u x can b r i n g about. The model c o n s i s t s of a system of f i v e r e a c t i o n - d i f f u s i o n equations d e s c r i b i n g the e v o l u t i o n of the v a r i o u s i o n i c c o n c e n t r a t i o n changes which are taken i n t o account. The i n f l u x of calcium ions through the nerve c e l l membrane r e s u l t s in n o n l i n e a r boundary c o n d i t i o n s which make a n a l y t i c a l s o l u t i o n u n obtainable. Through simple mathematical manipulations, we reduced the model to a system of three r e a c t i o n - d i f f u s i o n equations and we found that the system possesses some i n t e r e s t i n g mathematical p r o p e r t i e s such as a c o n s e r v a t i o n law and an i n v a r i a n t r e g i o n . The c o n s e r v a t i o n law p l a y s a key r o l e i n e s t a b l i s h i n g the uniqueness and asymptotic s t a b i l i t y of the s t e a d y - s t a t e which was shown to be s p a t i a l l y homogeneous and an a n a l y t i c a l e x p r e s s i o n for the s t e a d y - s t a t e value i s gi v e n . On the other hand, the i n v a r i a n t r e g i o n simply confirms the chemical law that the c o n c e n t r a t i o n l e v e l s of the complexes at a l l times must be gre a t e r than the i n i t i a l c o n c e n t r a t i o n l e v e l s but l e s s than the i n i t i a l c o n c e n t r a t i o n l e v e l s p l u s the t o t a l amount loaded i n t o 57 the system. Fu r t h e r examination of the equations r e v e a l s that the system possesses two small parameters, e and 5, whose presence a l l o w s us to c o n s t r u c t an approximation to the i n i t i a l t r a n s i e n t s o l u t i o n using a r e g u l a r p e r t u r b a t i o n technique. However, because of the n o n l i n e a r boundary c o n d i t i o n s , we were not able to f i n d an a n a l y t i c a l approximation to the i n i t i a l t r a n s i e n t s o l u t i o n which i s v a l i d f o r a longer p e r i o d of time. To t h i s end, a s i m p l i f i e d comparison theorem was used to provide rough bounds on the exact s o l u t i o n s . In order to give a complete d e s c r i p t i o n of the c a l c i u m dynamics, a B - s p l i n e c o l l o c a t i o n code, c a l l e d PDECOL, i s used to s o l v e the f u l l system of equations. The r e s u l t s are i d e n t i c a l to those given by Connor and Nikolakopoulou [13]. However we b e l i e v e that our numerical scheme i s more e f f i c i e n t than the POST a l g o r i t h m used in [13]. In c l o s i n g , we note that our a n a l y t i c a l r e s u l t s i n t h i s t h e s i s have pr o v i d e d g r e a t e r understanding and i n s i g h t i n t o the mechanisms of d i f f u s i o n , k i n e t i c s , and s e q u e s t r a t i o n of c a l c i u m in nerve c e l l s . 58 REFERENCES [I] Ahmed, Z., and Connor, J . A. (1979). "Measurement of c a l c i u m i n f l u x under v o l t a g e clamp in molluscan neurons using the metallochromic dye arsenazo I I I . " J . P h y s i o l (London) 286, pp. 61-82. [2] Ahmed, Z., and Connor, J . A. (1980). "Measurement of i n t e r n a l C a 2 + and r a p i d C a 2 + b u f f e r c a p a c i t y of molluscan neurons." Fed. Proc. V o l 39, page 2038. [3] Amundson, N.R. (1965) "Some f u r t h e r o b s e r v a t i o n s on t u b u l a r r e a c t o r s t a b i l i t y . " Can. J . chem. Engng 43, p. 49. [4] Aronson, D.G. and Weinberger, H.F. (1975) "Nonlinear d i f f u s i o n in p o p u l a t i o n g e n e t i c s , combustion and nerve pulse propagation", P a r t i a l D i f f e r e n t i a l Equations and R e l a t e d T o p i c s , L e c t u r e Notes i n Mathematics, S p r i n g e r -V e r l a g , New York, pp. 5-49. [5] Baker, P.F. Hodgkin, A.L., and Ridgway, E.B. (1971). " D e p o l a r i a z a t i o n and c a l c i u m e n t r y i n squid g i a n t axons." J . P h y s i o l . 218, pp. 709-755. [6] B l i n k s , J.R., Prendergast, F.G., and A l l e n , D.G. (1976). "Photoproteins as b i o l o g i c a l c a l c i u m i n d i c a t o r s . " Pharmacol. Rev. 28, pp. 1-93. [7] Boyce, W.E. and Diprima, R.C. (1977) Elementary  P i f f e r e n t i a l Equations and Boundary Value Problems. John Wiley and Sons, Inc. [8] Brown, J.E., Brown, P.K. and P i n t o , L.H. (1977). "Detection of l i g h t induced change of i n t r a c e l l u l a r i o n i z e d calcium c o n c e n t r a t i o n i n l i m u l u s v e n t r a l p h o t o r e c e p t o r s using arsenazo I I I . " J . P h y s i o l . (London) 267, pp. 299-320. [9] Brown, J.E., Cohen, L.B., DeWeer, P., P i n t o , L.H., Ross, W.N., and S a l z b e r g , B.M. (1975), "Rapid changes of i n t r a c e l l u l a r f r e e c a l c i u m c o n c e n t r a t i o n : d e t e c t i o n by metallochromic i n d i c a t o r dyes in squid g i a n t axon." Biophys. J . 15, pp. 1155-1160. [10] Chandra, J . and Davis, P.W. (1979) "Comparison theorems for systems of r e a c t i o n - d i f f u s i o n equations". A p p l i e d Nonlinear  A n a l y s i s V. Lakshmikantham (ed.), Acadamic Pr e s s , 1979, New York, N.Y. pp.79-88. [ I I ] Chueh, K., Conley, C , and Smoller, J . (1977) " P o s i t i v e l y i n v a r i a n t regions f o r systems of n o n l i n e a r equations". Indiana Univ. Math. J . , 26, pp. 373-392. 59 [12] Cohen, D.S., " M u l t i p l e s o l u t i o n s of n o n l i n e a r p a r t i a l d i f f e r e n t i a l equations", Nonlinear Problems i n the P h y s i c a l S c i e n c e s and B i o l o g y , L e c t u r e Notes in Mathematics, S p r i n g e r - V e r l a g , New York, 1972. [13] Connor, J.A. and Nikolakopoulou, G. "Calcium D i f f u s i o n and B u f f e r i n g i n Nerve Cytoplasm". L e c t u r e s on Mathematics in the L i f e S c i e n c e s , v o l . 15 (1982), pp. 79-102. a [14] Conway, E., Hoff, D., and Smoller, J . (1978) "Large time behaviour of s o l u t i o n s of systems of n o n l i n e a r r e a c t i o n -d i f f u s i o n equations". SIAM J . Appl. Math. 35 pp. 1-15. [15] Crank, J . (1975) The Mathematics of D i f f u s i o n . Clarendon p r e s s , Oxford. [16] E b a s h i , S. and Endo, M. (1967) "Calcium ion and muscle c o n t r a c t i o n " . Prog. Biophys. molec. B i o l . V o l 18, pp. 123-183. [17] Jackson, R.(1973) "A simple geometric c o n d i t i o n f o r i n s t a b i l i t y i n c a t a l y s t p e l l e t s at u n i t Lewis number." Chem. Engng S c i . 28, p. 1355. [18] J o b s i s , F.F. and O'Connor, M.J. (1966) "Calcium r e l e a s e and r e a b s o r p t i o n i n the s a r t o r i u s muscle of the toad." Biochem. biophys. Res. Commun. 25, pp. 246-252. [19] K r a s n o s e l ' s k i i , M. A. (1964) T o p o l o g i c a l methods in the  theory of n o n l i n e a r i n t e r n a l equations. Pergamon press, Oxford. [20] Madsen, N.K. amd Somcovec, R.F. (1976) "General Software f o r P a r t i a l D i f f e r e n t i a l E q u a t i o n s " . Numerical Methods for  D i f f e r e n t i a l Systems j_ Recent Developments i n Algorithms,  Software, and A p p l i c a t i o n s , L. Lapidus and W.E. S c h i e s s e r , e d i t o r s , Academic Press Inc., New York, pp. 229-242. [21] M i l e s , F.A. Exc i t a b l e C e l l s . (1969) W i l l i a m Heinemann Med i c a l Books L i m i t e d , London. [22] Moore, C. (1982) "The s o l u t i o n of P a r t i a l D i f f e r e n t i a l Equations i n one s p a t i a l dimension". UBC PDECOL documentation. [23] Sandow, A. (1965) " E x - c i t a t i o n - c o n t r a c t i o n c o u p l i n g in s k e l e t a l muscle." Pharmac. Rev. 17, pp. 265-320. [24] Scarpa, A., B r i n l e y , F . J . , T i f f e r t , T. and Dubyak, G.R. (1978). "Metallochromic i n d i c a t o r s at i o n i z e d c a l cium." Ann. N^ Y^ Acad. S c i 307 , pp. 86-112. [25] Selby S.M. (1967) Standard Mathematical Tables. The Chemical Rubber Company. 60 APPENDIX A : DISCUSSION OF THE NUMERICAL ALGORITHM - PDECOL The UBC l i b r a r y program, PDECOL, i s used to solve the f u l l system of equations (3.2)-(3.4). PDECOL i s a general purpose program f o r n u m e r i c a l l y s o l v i n g systems of p a r t i a l d i f f e r e n t i a l e quations i n one s p a t i a l dimension. The c l a s s of problems which PDECOL can solve must have the f o l l o w i n g s t r u c t u r e : 3u/9t = F ( x , t , u , 9 u / 9 x , 9 2 u / 9 x 2 ) , (A.1) where u = (u ,...,u ), F = (f , . . . , f ), x e I = [a,b] and t > t 0 . 1 n i n The boundary c o n d i t i o n s must have the form B(u,3u/3x) = z ( t ) , x e 31, t > t 0 , (A.2) where B = (b ,...,b ) and z = (z ,...,z ). The i n i t i a l c o n d i t i o n s 1 n 1 n on a ^ x < b are u ( x , t 0 ) = g ( x ) , (A.3) where g = (g ,...,g ). 1 n The i n i t i a l c o n d i t i o n s must be c o n s i s t e n t with the boundary c o n d i t i o n s , i . e . , B(g,3g/3x) = z ( t 0 ) , x e 31. (A.4) F i n a l l y a l l f u n c t i o n s must be continuous i n t and piecewise continuous i n x. Given a system of PDEs of the above s t r u c t u r e , PDECOL s o l v e s the system using time d i s c r e t i z a t i o n and c o l l o c a t i o n s p a t i a l 61 d i s c r e t i z a t i o n t e c h n i q u e s 1 . For s i m p l i c i t y , we s h a l l set n = 1 and o u t l i n e the s o l u t i o n procedures used by PDECOL. For more d e t a i l e d i n f o r m a t i o n we r e f e r the reader to [22]. PDECOL Al g o r i t h m The a l g o r i t h m r e q u i r e s the user to s p e c i f y a s p a t i a l mesh (x ,...,x ) on [a,b] such that 1 N a = x < x < . . . < x =b . 1 2 N A s s o c i a t e d with each i n t e r v a l i n the mesh are the polynomials p ( x ) , i = 1,...,N-1, which form a piecewise polynomial space i t h at i s used to compute the approximate s o l u t i o n . The user must s p e c i f y k, the order of the polynomials (k = the polynomial degree + 1) and the number of c o n t i n u i t y c o n d i t i o n s (NCC) to be imposed on the polynomial p i e c e s a c r o s s the mesh p o i n t s x . The i recommended c h o i c e s are k = 6 and NCC = 2 so that the approximate s o l u t i o n which i s made up of separate f i f t h - d e g r e e polynomial p i e c e s and the f i r s t d e r i v a t i v e of the approximate s o l u t i o n are both continuous across the mesh p o i n t s and hence on the e n t i r e i n t e r v a l [ a , b ] . The dimension of the l i n e a r polynomial space ( i . e . , the t o t a l number of degrees of freedom) over [a,b] i s t h e r e f o r e m = k(N-1) - NCC(N-2). 1 T h e o r e t i c a l r e s u l t s [20] suggest that the p r o v i d e s more r e l i a b l e r e s u l t s than f i n i t e - d i f f e r e n c e method. c o l l o c a t i o n method a s t r a i g h t f o r w a r d 62 In PDECOL a B - s p l i n e b a s i s c o n s i s t i n g of m known piecewise polynomial f u n c t i o n s 6 (x) are used to span the polynomial space. i T h i s b a s i s has the pr o p e r t y that f o r any x e [ a , b ] , at most k of the 6 (x) have nonzero v a l u e s . T h i s property r e s u l t s i n a banded i m atrix which leads to e f f i c i e n t numerical computation. For a given value of t , the a l g o r i t h m assumes an approximate s o l u t i o n of the form m U(x,t) = Z c ( t ) 0 (x) , i = 1 i i where the c o e f f i c i e n t s c ( t ) , i = 1,...,m, are determined by a i s p a t i a l c o l l o c a t i o n technique. The c o l l o c a t i o n technique r e q u i r e s U(x,t) to s a t i s f y the system of PDEs and the boundary c o n d i t i o n s at a set of m c o l l o c a t i o n p o i n t s . The c o l l o c a t i o n p o i n t s , £ , are j s e l e c t e d a u t o m a t i c a l l y by PDECOL such that a = £ < £ < . . . < £ = b . 1 2 m Hence f o r each j = 1,...,m, m u(* ,t) = U U ,t) = Z c (t)0 U ), (A.5) j j i=1 i i j where u i s the exact s o l u t i o n to the PDE. S u b s t i t u t i n g equation (A.5) i n t o (A.1) we o b t a i n a system of time dependent o r d i n a r y d i f f e r e n t i a l equations m Z 6 (£ )9c ( t ) / 9 t = F(£ ,t,UU ,t),9U(£ ,t)/3x,9 2uU , t ) / 3 x 2 ) , i=1 i j i j j j j (A.6) 6 3 s u b j e c t to the i n i t i a l c o n d i t i o n s m I c ( t o ) 0 U ) = g U ) . (A.7) i=1 i i j j Equations (A.6)-(A.7) must be s a t i s f i e d f o r j = 2,...,m-1. At the boundaries, s p e c i a l c o l l o c a t i o n equations are formed to account for the boundary c o n d i t i o n s . Thus we have transformed a s i n g l e PDE i n t o a system of o r d i n a r y d i f f e r e n t i a l equations which can be solved by i n t e g r a t i n g with respect to the time v a r i a b l e . The present a l g o r i t h m possesses three d e s i r a b l e f e a t u r e s . F i r s t of a l l , there are no l i m i t s on the number of PDEs to be so l v e d . Secondly, the s o l u t i o n procedures a p p l i e d to the system of time dependent o r d i n a r y d i f f e r e n t i a l equations (A.6)-(A.7) are very e f f i c i e n t s i n c e the cho i c e of the B - s p l i n e b a s i s r e s u l t s i n a banded mxm matrix of maximum bandwidth [ 2 n ( k - 1 ) - l ] , T h i r d l y , PDECOL allows the user to r e s t a r t i n t e g r a t i o n from the prev i o u s t e r m i n a t i o n time rather than to repeat the e n t i r e run. T h i s i s p a r t i c u l a r l y convenient and economical s i n c e the length of time to reach the e q u i l i b r i u m s t a t e i s not known a p r i o r i . D e spite a l l of the above d e s i r a b l e f e a t u r e s , PDECOL s u f f e r s from one major drawback 1, namely, the c o n s i s t e n c y c o n d i t i o n (A.4) must be s a t i s f i e d by the PDEs i n order to invoke PDECOL c o r r e c t l y . T h i s i s a severe r e s t r i c t i o n s i n c e i n a p p l i c a t i o n s 'Through p r i v a t e c o n s e r v a t i o n s with Carolyn Moore and Tom N i c o l , both of whom are numerical a n a l y s t s at the UBC Computer Center, we are unaware of any other c o l l o c a t i o n method which w i l l a u t o m a t i c a l l y handle the i n c o n s i s t e n c y s i t u a t i o n . 64 where an impulse or an i n f l u x i s in t r o d u c e d i n t o the system at the boundaries at t = t 0 + , the c o n s i s t e n c y requirement on the PDEs i s o f t e n v i o l a t e d . In our case the i n i t i a l and boundary c o n d i t i o n s on V(x,t) are V(x,0) = V 0, r A + B l n V ( 0 , t ) , 0 < t < t * , 9V(0,t)/3x = L 0, t > t * . Now s i n c e l n V 0 * -A/B, 3V(x,0)/3x * 3V(0,t)/9x at (x,t) = (0,0), and thus the c o n s i s t e n c y c o n d i t i o n i s not s a t i s f i e d . Moreover, we observe that 9V(0,t)/3x i s not continuous at t = t * . In order to overcome our present numerical d i f f i c u l t y , we intr o d u c e d a t r a p e z o i d a l - s h a p e d f u n c t i o n H(t) = 1/2 [ t a n h ( 7 ( t - t ' ) ) - t a n h ( 7 ( t - t " ) ) ] where we have chosen 7 = 1 0 1 6 , t' = 10" 1 ", t " = t * - t ' . The above parameter values are chosen such that H(0) * 0 ~ H ( t * ) , H(t') ~ 1/2 * H ( t " ) , H(t) - 1, 2t' < t < t * - 2 t ' , where the "-" sign i n d i c a t e s an e r r o r of l e s s than 10~*°. A sketch of H(t) i s shown i n F i g u r e 6. 65 H ( t ) 1.0-0 . 5 -0.0 2 t ' t * - 2 t ' t * F i g u r e 6 : Sketch of the t r a p e z o i d a l - s h a p e d f u n c t i o n H ( t ) . To overcome our numerical d i f f i c u l t y , we can reformulate our boundary c o n d i t i o n s as Now s i n c e H(0) = 0 = H ( t * ) , the c o n s i s t e n c y c o n d i t i o n i s s a t i s f i e d and 3V(0,t)/9x i s continuous everywhere. Although the above method i s formulated to s u i t our p a r t i c u l a r s i t u a t i o n , i t i s c l e a r t h a t the method can be m o d i f i e d and adapted to other s i m i l a r s i t u a t i o n s as w e l l . We s h a l l devote the remainder of t h i s Appendix to d i s c u s s the numerical accuracy of PDECOL. There are mainly two sources of e r r o r i n the approximate s o l u t i o n generated by PDECOL. The f i r s t i s due to the time d i s c r e t i z a t i o n method and the second i s due to B l n V ( 0 , t ) ] H ( t ) , 0 < t < t * , t > t * . 66 the c o l l o c a t i o n s p a t i a l d i s c r e t i z a t i o n technique. PDECOL c o n t r o l s the time d i s c r e t i z a t i o n e r r o r by s e l e c t i n g a s t e p s i z e below the user's s p e c i f i c a t i o n . However, i t has no c o n t r o l over the s p a t i a l d i s c r e t i z a t i o n e r r o r s . T h e o r e t i c a l r e s u l t s i n [20] i n d i c a t e that the maximum e r r o r between the exact s o l u t i o n and the approximate s o l u t i o n generated k by PDECOL i s p r o p o r t i o n a l to h , where h = max |x - x |. i i i-1 Moreover the maximum e r r o r 1 between the j t h d e r i v a t i v e s of the exact s o l u t i o n and that of the approximate s o l u t i o n i s k-j p r o p o r t i o n a l to h . Our system of r e a c t i o n - d i f f u s i o n equations i n v o l v e s second order s p a t i a l d e r i v a t i v e s and hence the maximum k-2 e r r o r i s of the order h 1 I n p r a c t i c e one should choose a v a r i a b l e s p a t i a l mesh so that more mesh p o i n t s are i n s e r t e d near the region where the s o l u t i o n changes r a p i d l y . 67 Appendix B : AN APPROXIMATION TO A SECOND-ORDER PDE In t h i s appendix we s h a l l f i n d an approximate s o l u t i o n to the l i n e a r second-order p a r t i a l d i f f e r e n t i a l equation dz/dt = a 9 2 z / 9 x 2 + bw(x,t) (B.I) with i n i t i a l c o n d i t i o n on 0 < x < L z(x,0) = 0, (B.2) and boundary c o n d i t i o n s 9z(0,t)/9x - a ( t ) z ( 0 , t ) = 0, 9z(L,t)/9x =• 0, (B .3) and where a ( t ) and w(x,t) are known f u n c t i o n s and a, b are given c o n s t a n t s . We s h a l l t e m p o r a r i l y set a ( t ) = a = constant and seek a s o l u t i o n to (B.1)-(B.3) of the form OO z ( x , t ) = I T (t)cosX ( x - L ) . (B.4) n= 1 n n In order to s a t i s f y the boundary c o n d i t i o n s , X must s a t i s f y n the t r a n s c e n d e n t a l equation X sinX L - acosX L = 0, (B.5) n n n for n = 1,2,... with 0 < X, < X 2 < X 3 < .... We expand OO w(x,t) = Z w (t)cosX (x-L) (B.6) n= 1 n n 68 where w (t) = n {/ w(x,t)cosX (x-L)dx}/K , X * 0, 0 n n n [J w(x,t)dx]/L, 0 X = 0, n f o r n = 1,2 , and where K = [L + acos 2X L/X 2]/2. ( Note : n n n X, = 0 i f f a = 0.) S u b s t i t u t i n g equations (B.4) and (B.6) i n t o (B.1)-(B.2) and s o l v i n g y i e l d s T (t) = b J w ( s ) e x p ( - a X 2 ( t - s ) ) d s . n O n n (B.7) Thus f o r con s t a n t a, a s o l u t i o n to the p a r t i a l d i f f e r e n t i a l equation (B.1)-(B.3) i s given by (B.4) with T (t) d e f i n e d by n (B.7) . The f o l l o w i n g lemma extends the above r e s u l t f o r general a ( t ) . Lemma 1 : Let a ( t ) > 0 be continuous and d i f f e r e n t i a b l e . Provided t h a t | A a / o | = |[a(t+At) - a ( t ) ] / a ( t ) | « 1, then (B.4) i s an a r b i t r a r i l y good approximation to the exact s o l u t i o n of the p a r t i a l d i f f e r e n t i a l equation (B.1)-(B.3). Proof : I t s u f f i c e s to show that |AX /X | = |[X (t+At) - X ( t ) ] / X ( t ) | < |Aa/a|, n n n n n for a l l n > 1 so that the roots {X } to the t r a n s c e n d e n t a l n 69 equation are approximately constant and hence equations (B .4)-(B.7) h o l d . We set a = a i t ) , X = X (t) and d i f f e r e n t i a t e (B.5) with n n respect to t g i v i n g X '(t)/X = {1/G(X ,a)}a'it)/a, n n n or |X'(t)/X | = |{1/G(X , a ) } a ' ( t ) / a | n n n where G(X ,a) = 1 + aL + X 2L/a. n n Now min |G(X ,a)| = |G(X ,X )| = |1 + 2X L|. a n n n n Thus |AX /X | < [1/(1 + 2X L)]|Aa/a| < |Aa/a| n n n for a l l n > 1 as d e s i r e d . The key step i n s o l v i n g the p a r t i a l d i f f e r e n t i a l equation (B.1)-(B.3) l i e s in f i n d i n g the roo t s of the t r a n s c e n d e n t a l equation (B.5). Numerical schemes w i l l be more e f f i c i e n t i f we can g i v e an approximate l o c a t i o n of the r o o t s . The next lemma serves t h i s purpose. Lemma 2 : Let j3 = nff/L, n = 0 , 1 , 2 , . . . n denote the roots to (B.5) with a = 0. For a small p o s i t i v e p e r t u r b a t i o n Aa, the roo t s to (B.5) are approximately given by X, * • (Aa/L) , X * 0 + Aa/( (n-1 )ir) , n > 1. n n-1 70 Proof : Let $ = X L and r e p l a c e a by Aa i n ( B . 5 ) . I t i s easy to n n deduce from equation (B .5) that tan$ = tan($ - (n-1 )ir) = AaL/S , n n n and thus $ - (n-1)7r = t a n " 1 ( A a L / $ ) « (AaL/$ ) n n n for |AaL/5 | << 1. The l a s t approximation i s j u s t i f i e d s i n c e the n p e r t u r b a t i o n Aa i s assumed to be small and |AaL/<; | decreases as n n i n c r e a s e s . Thus for n = 1, $, ~ (AaL/$ , ) , or X, * v/(Aa/L),. and f o r n > 1, 5 ~ . ( n - l ) i r + (AaL/$ ) or X = 0 + A a / ( ( n - l ) i r ) n n n n-1 as r e q u i r e d . 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0079345/manifest

Comment

Related Items