META-STABILITY OF T H E GIERER MEINHARDT EQUATIONS by DAVID IRON B.A.Sc. (Mechanical Engineering) University of Toronto, 1988 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF M A S T E R OF SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Mathematics Institute of Applied Mathematics We accept this thesis as conforming to the requireqVstandard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A September 1997 © David Iron, 1997 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer-ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada A b s t r a c t A well-known system of partial differential equations, known as the Gierer Meinhardt system, has been used to model cellular differentiation and morphogenesis. The system is of reaction-diffusion type and involves the determination of an activator and an in-hibitor concentration field. Long-lived isolated spike solutions for the activator model the localized concentration profile that is responsible for cellular differentiation. In a biological context, the Gierer Meinhardt system has been used to model such events as head determination in the hydra and heart formation in axolotl. This thesis involves a careful numerical and asymptotic analysis of this system in one dimension for a specific parameter set and a limited analysis of this system in a multi-dimensional setting. Numerical analysis has revealed that once the spikes form they continue to move on an extremely slow time scale. This type of phenomenon is a general indicator of meta-stable behaviour. By perturbing off of an isolated spike solution an exponentially small eigenvalue of the linearized operator was found. This small eigenvalue accounted for the extremely slow motion found numerically and thus was used to obtain an equation of motion for the location of the spike. The Gierer Meinhardt system is analyzed in the limit of small activator diffusivity for both a finite inhibitor diffusivity and for an asymptotically large inhibitor diffusivity. In this thesis, the mathematical techniques used include the method of matched asymptotic expansions, spectral theory and numerical computations. i i T a b l e o f C o n t e n t s Abstract ii Table of Contents iii Acknowledgement iv Chapter 1. Introduction 1 Chapter 2. Infinite Inhibitor Diffusion Coefficient 12 2.1 Introduction 12 2.2 A One-Spike Quasi-Equilibrium Solution 14 2.3 The One-Spike Linear Eigenvalue Problem 16 2.4 A n Exponentially Small Eigenvalue 24 2.5 The Slow Motion of the Spike 29 2.6 A n n-Spike Solution 31 Chapter 3. Finite Inhibitor Diffusion Coefficient 38 3.1 A One-Spike Quasi-Equilibrium Solution 38 3.2 A One-Spike Eigenvalue Problem 45 3.3 A n n-Spike Solution 49 3.4 The n-Spike Linearized Eigenvalue Problem 54 Chapter 4. A Spike in a Multi-Dimensional Domain 62 4.1 A n Exponentially Small Eigenvalue 68 Chapter 5. Conclusions 78 Bibliography 81 i i i A c k n o w l e d g e m e n t This thesis was produced under the supervision of Dr. Michael Ward and Dr. Robert Miura. I would like to thank them both for presenting the problem to me as well as providing advice, encouragement and support during the preparation of this thesis. iv C h a p t e r 1 I n t r o d u c t i o n The development of a complete organism from a single cell is still one of the great mys-teries remaining in biological science. Many different mechanisms are involved in the completion of this process. Some involve mechanical interactions between the cells, and between the cells and their extracellular matrix, such as gastrulation. In other process such as organogenesis, among a group of similar cells, certain cells will become differen-tiated from their neighbors. These cells will begin to change and develop the necessary structures for the organs that they will eventually form. The mechanism responsible for cell differentiation varies for different structures. Experiments have shown that a local increase in the concentration of a substance called a morphogen, or inducer, is often responsible for organogenesis. The inducer will cause the activation of genes which will then produce the specific proteins used by the mature organ. Thus, cells in the neigh-borhood of an inducer concentration peak will form one organ and the surrounding cells will have other fates. In some cases, isolates spikes are required, as in the formation of the heart or liver. In other cases, such as the spinal cord, the periodic nature of the resulting structure would require periodic fluctuations of the activator. In all cases, precise positioning of the structure is required for the resulting organism to be viable. The mechanism for placement of the concentration spike must be stable to the random fluctuations present in any biological system. Turing [13] proposed a reaction-diffusion system of activator-inhibitor type that suggested 1 Chapter 1. Introduction that a two species chemical system with Fickian diffusion and non-linear reactive terms could model morphogenesis. He conjectured that some stable spatially inhomogeneous solutions to this system could have isolated peaks in the inducer concentration. As a first step to exploring this hypothesis, he examined the spectrum of the linearized reaction-diffusion system about a spatially homogeneous equilibrium solution. He found that, under certain constraints, a finite number of spatially periodic eigenmodes will have positive eigenvalues. Subsequent studies (e.g. Gierer and Meinhardt [5], Holloway [6]), which have involved large-scale numerical computations, have shown that these eigenmodes will grow in time until they enter the non-linear regime. Nonlinear effects will then lead to a saturation of the amplitudes of these modes. When this occurs, isolated spikes of the activator concentration will typically be formed. A qualitative explanation for this phenomenon is as follows. The activator is auto-catalytic, and the inhibitor diffuses rapidly and slows the production of the activator. It itself is catalyzed by the activator. Any local increase in the activator concentration will continue to increase due to auto-catalysis. This, eventually, will lead to the formation of a spike. The local increase in the activator concentration will cause a local increase in the inhibitor concentration, which will then spread quickly. This globally elevated concentration of the inhibitor will localize the existing spike and will also prevent the formation of additional spikes in the activator concentration at other spatial locations. In this thesis we analyze spike behavior for the following general Gierer Meinhardt system in one spatial dimension. In this system, the activator concentration A = A(x, t) and the inhibitor concentration H = H(x,t) satisfy Ap At = DaAxx - ^ A + paCa—, -L<x<L, t>0, (1.1a) Am Ht = DhHxx- \xhE + Chph—, -L<x<L, t>0, (1.1b) Hx(±L,t) = 0, Ax(±L,t)=0. (1.1c) 2 Chapter 1. Introduction The exponents p, q, m, and s are assumed to satisfy v — 1 ui p > 1, q > 0, m > 0, s > 0, 0 < < (1-2) 9 5 + 1 V ' The values of (p, q, m, s) will depend on the details of the reaction. The constant pa represents the rate of increase in active sources caused by the presence of activator and inhibited by the presence of inhibitor. The constant ph is the rate of increase of active sources of the inhibitors that are turned on by the activator. In addition, Da and are the diffusion coefficients of the activator and the inhibitor, respectively, and and Ca and Ch are the coupling constants. The parameter set (p,q,m,s) — (2,1,2,0) is used to model a system in which the activator and inhibitor have different sources. The set (p, q, m, s) = (2, 4, 2,4) is used to represent an activator-inhibitor system with common sources. Gierer and Meinhardt proceeded to use these equations to model the head formation in the hydra. We may reduce the number of parameters appearing in the Gierer Meinhardt system using an appropriate non-dimensionalization of the problem. We choose, t = t'T, A = A0A, H = H0H, x = Lx', (1.3) where, PhChA™ A0 X = ChPh PaCg s+1 Y 1 f PaPa ChPh) Pa P-a qm- (p- l)(s + 1)' This results in the following non-dimensional system: H ^(±1 ,0 = 0, ^(±i,t') = o. Av = D'aAx,x, -A + rhHf — D'hHxixi — pH + Ph rn -1<X'<1, t'>0, -l<x'<l, t'>0, (1.4) (1.5) (1.6) (1.7a) (1.7b) (1.7c) 3 Chapter 1. Introduction Here Da Dh p = phH0. (1.8) L2Pa Bifurcation and perturbation methods have been the two main analytical methods used to examine the behavior of solutions to the non-linear Gierer Meinhardt system. Previous analyses have shown that when the diffusion coefficients are sufficiently large, the spa-tially homogeneous solution is stable. As one or both of the diffusion coefficients become smaller, this solution becomes unstable and spike-like patterns in the activator concen-tration will result. Bifurcation analysis is used to investigate the properties of solutions near this bifurcation point. In general, the limitation of this method is that it will lead only to small amplitude solutions that bifurcate off of the trivial solution. However, it is the large amplitude solutions which are of interest in morphogenesis. The calculation of these solutions typically requires a full numerical simulation. In certain cases, pertur-bation methods have been used to calculate large amplitude solutions. Keener[7] used perturbation methods to investigate the nature of large amplitude steady-state spike so-lutions in the limit for which the diffusion coefficient of the inhibitor tends to infinity. This analysis leads to the non-local problem studied in the second chapter of this thesis. The analysis done by Nishiura[10] links the bifurcation analysis and the perturbation analysis. Before we describe the goals and the outline of the thesis and summarize some previous work, we find an appropriate scaling of (1.7) for spike solutions. We introduce a small parameter e in (1.8) by In the variables of (1.7) the amplitude of a spike solution tends to infinity as e —> 0. Therefore, it is convenient to introduce new variables so that the amplitude of the spike 4 Chapter 1. Introduction solution is 0(1) as e —>• 0. For simplicity, in what follows, we drop the primes in (1.7). We first introduce a and h by A = e~Uaa, H = e-"hh, (1.10) where the exponents va and vh are to be found. To balance the terms in (1.7a) we require, -va =-VaP + qvh- (1-11) We are interested in solutions involving isolated spikes of the activator concentration. We therefore expect A to be localized to within an 0(e) region near the spike. Thus in our scaling of (1.7b) we will consider an averaged balancing. Specifically, we integrate (1.7b) over the domain to get /l rl rl Am Htdx = -pJ Hdx + J -jjjdx. (1.12) Since A will be localized to within an 0(e) region about the spike location x 0 , we scale x in the last term by y = (x — x 0 ) e _ 1 . Balancing the terms in this equation results in the following: -vh = -uam + uhs + 1. (1.13) Solving equations (1.11) and (1.13) yields, (p - l)(s + 1) — qm' (p - l)(s + 1) — qm ' This determines the scaling in (1.10). In terms of these new variables, (1.7) becomes a? at = e2axx — a + — , —1 < x < 1, t>0, (1.15a) am rhht = Dhhxx - ph + e'1-—, - 1 < x < 1, t>0, (1.15b) ns a x ( ± l , t ) = 0 , hx(±l,t) = 0. (1.15c) 5 Chapter 1. Introduction We will study this scaled system analytically and numerically as e -» 0 for two different ranges of Dh: Dh —> oo , weak coupling limit; Dh = 0(1) , strong coupling limit. (1-16) In the weak coupling limit, the inhibitor will diffuse very rapidly compared to the size of the domain. Thus, the concentration of the inhibitor may be considered to be constant in space. Each spike in the activator concentration is confined to a small region and will thus act as a point source of inhibitor. The equilibrium level of inhibitor concentration will then, in effect, count the number of spikes of activator concentration and the position of the spikes will be irrelevant. Too many spikes will cause the equilibrium level of inhibitor to become large and the spikes will become unstable. In the strong coupling limit, the inhibitor will still diffuse much faster then the activator, but the length scale of its diffusion is comparable to the size of the domain. Thus, if the distance between adjacent spikes is small, large levels of inhibitor concentration may build up in this area. However, when the distance between adjacent spikes is large, the spikes do not feel each others presence, since the inhibitor concentration decays exponentially with the distance from the source of inhibitor. Therefore, in this strong coupling limit, the positioning of the spikes will play an important role in determining the stability of a configuration of spikes. Previous work on the Gierer Meinhardt system has focused on small amplitude solu-tions. In this thesis, we will attempt to construct large amplitude equilibrium and time-dependent solutions. The analysis will be done for the limit e —>• 0 for two different ranges of Dh (Dh —>• oo in chapter 2 and Dh = 0(1) in chapter 3). Our preliminary numerical computations have suggested that spike solutions to the Gierer Meinhardt system will be formed quickly in time from initial data. These spike solutions persist in their basic shape, but the centers of the spike layers migrate very slowly towards their equilibrium positions. This type of phenomenon, in which internal layers move exceedingly slowly in 6 Chapter 1. Introduction time, is referred to as meta-stable behavior. The motivation of this study is the numerical simulations presented in the thesis of David M . Holloway[6]. The parameter values used in this thesis correspond to the strong coupling limit where Dh = 0(1). In Holloway's thesis, numerical simulations using a finite difference method were run from 20 000 to 560 000 iterations at a fixed time step before an equilibrium was achieved for the discretized problem. This very slow convergence of the system towards equilibrium suggests that the system could exhibit meta-stable behavior. Simulations carried out in two dimensions resulted in a somewhat random pattern of equilibrium spike positions in the computed solution. I believe that the randomness of the spike locations for the computed equilibrium solutions does not correspond to a true equilibrium solution for (1.7), but is instead likely due to meta-stable behavior of some quasi-equilibrium solution. Since meta-stable solutions evolve on such a slow time scale, these quasi-equilibrium solutions could easily be mistaken for true equilibrium solutions. In a one dimensional domain, true equilibrium solutions have equally spaced spike locations. It is conjectured that the analogous result, in a two dimensional domain, is that an equilibrium spike layer solution should have spikes that lie on lattice sites and not on random positions in the domain. Our goal is to ascertain if meta-stable behavior occurs for (1.7). Meta-stability has been studied previously for other partial differential equations (e.g. Ward [17]). As shown in this previous work, a necessary condition for meta-stability is that the spectrum of the linearization of the partial differential equation about some canonical spike-type or shock-type profile contains asymptotically exponentially small eigenvalues in the limit for which the width of the spike or shock profile tends to zero. The existence of these eigenvalues is usually indicated by a near indeterminacy in determining internal layer locations corresponding to certain equilibrium solutions. 7 Chapter 1. Introduction To illustrate this phenomena consider the following two problems on |x| < 1, t > 0: ut = e2uxx + 2(u-u3), u s ( ± l ) = 0, (1.17) ut = e2uxx - u + u2, ux(±l) = 0. (1-18) Equation (1.17) is a phase transition problem, which gives rise to shock solutions. Equa-tion (1.18) resembles the activator equation when the inhibitor is a given constant. The canonical one-shock profile for (1.17) has the form us(y) = tanh(y). Consider the function UE{X) = us (£=£a) that satisfies the steady-state equation corresponding to (1.17). Here XQ is a constant satisfying |x 0 | < 1- Since this function fails to satisfy the boundary condition in (1.17) by only exponentially small terms for any rr0 in \XQ\ < 1, it is analytically very difficult to determine the correct value x0 = 0 corresponding to a true equilibrium solution. Hence, we shall refer to UE(X), where XQ is arbitrary in |x 0 | < 1, as a quasi-equilibrium solution. To link this near indeterminacy to the occurrence of an exponentially small eigenvalue, we linearize (1.17) about our quasi-equilibrium solution uE(x). This leads to the eigenvalue problem L(j) = t2(f)xx + (2 - 6u2E)(/) = A<6, - 1 < x < 1, (1.19) </>x(±l) = 0. (1.20) Since UE solves the steady problem for (1.17) it follows that Lu'E = 0. Hence, on the infinite line —oo < x < oo subject to (j) —>• 0 as x —>• ±oo we have that <f> = u'E and A = 0 is an eigenpair of (1.19). However, since u'E fails to satisfy the boundary conditions for the finite domain problem by only exponentially small amounts, we expect that the finite boundary will perturb this eigenpair by only exponentially small terms. Therefore, this suggest that there is an eigenvalue of (1.19) which is exponentially small. Moreover, since u's > 0, it follows that u'E has no nodal points. Hence the exponentially small eigenvalue must be the principal eigenvalue. It is this eigenvalue that is responsible for the meta-stable behavior that occurs for the corresponding time-dependent problem. As 8 Chapter 1. Introduction a remark, a similar situation arises for a solution with n shock layers. In this case, the quasi-equilibrium solution has the form unE(x) = Y^=ius {s~T^)i f ° r some X{ satisfying \xi\ < 1. The eigenvalue problem associated with the linearization of (1.17) about unE has n exponentially small eigenvalues, one associated with each internal layer. These n exponentially small eigenvalues lead to the slow coupling between shock layers for the evolution problem. For a precise quantitative description of these results see the references in [17]. A similar analysis may be applied to (1.18). Here the canonical spike profile is given by us(x) = |sech 2 ( | ) . Again the quasi-equilibrium solution uE = us ( 5-^ f l) will satisfy the steady-state equation corresponding to (1.18) but fails to satisfy the boundary conditions in (1.18) by only exponentially small amounts for any value of XQ in \XQ\ < 1. Thus, determining the true equilibrium value xo = 0 requires exponential precision. Linearizing (1.18) about uE results in the eigenvalue problem, L<j> = e2(bxx + (-1 + 2uE)<f> = \<f>, (1.21) 0x(±l) = 0. (1.22) It is clear that Lu'E = 0 and that u'E fails to satisfy the Neumann boundary conditions in this problem by only exponentially small amounts. Thus, there must be an eigenpair exponentially close to A = 0 and (j) = u'E. This case differs from the shock problem (1.17) in that now u'E has exactly one nodal point. Therefore, uE must be exponentially close to the second eigenfunction of (1.21). Thus, the exponentially small eigenvalue is not the principal eigenvalue for (1.21) and hence there is no reason to expect that meta-stability will occur for (1.18). This suggests that the Gierer Meinhardt equations, under the as-sumption that h is a given constant, may not exhibit meta-stable behavior. We will show that meta-stable behavior results from the coupling of the activator and inhibitor con-centration fields. We will also show that there are exponentially small eigenvalues for the activator-inhibitor problem and that, under appropriate conditions, these exponentially 9 Chapter 1. Introduction small eigenvalues are indeed the principal eigenvalues. There are also a few rigorous results for the Gierer Meinhardt in certain limiting situa-tions. Multi-peak equilibrium solutions to the Gierer Meinhardt equations are rigorously shown to exist in one-dimensional domains [12]. Similar results for multi-dimensional domains can be found in [9]. These papers provide interesting examples of rigorous existence results as they also provide a qualitative description of the solutions. The organization of this thesis is as follows. In Chapter 2 we will consider the weak coupling limit Dh —>• oo. This leads to the what is known as the Shadow system introduced in [10]. A one-spike quasi-equilibrium solution to the Shadow system will be constructed using the method of matched asymptotic expansions. The eigenvalue problem associated with the linearization about this solution will be obtained. The spectrum of this problem will then be examined and an exponentially small eigenvalue will be shown to exist. Under some appropriate conditions, this eigenvalue will be demonstrated to be the principal eigenvalue. Then, the analysis of metastable behavior associated with phase transition problems considered in [17] will be extended to quantify the meta-stable behavior in our system. This analysis, which is based on the projection method of [17], imposes a limiting solvability condition to derive an ordinary differential equation governing the motion of the center of one spike. Multiple spike solutions will then be considered. A similar spectral analysis to that of the one spike case, will reveal that the principal eigenvalue will not be exponentially small. Thus, solutions with multiple spikes are not meta-stable. In Chapter 3, we will consider the strong coupling case for which the inhibitor diffusion coefficient Dh is 0(1). The study of this case is significantly more intricate than the previous case in that we no longer have the simplified Shadow system to work with. Again we will use the method of matched asymptotic expansions to construct a one-spike quasi-equilibrium solution. In this case, the inhibitor concentration is no longer 10 Chapter 1. Introduction spatially constant. We will use the one-spike quasi-equilibrium solution to derive an eigenvalue problem as in the second chapter. This eigenvalue problem will prove to be of a similar form to the eigenvalue problem of the second chapter and thus the previous results may be applied. The n-spike quasi-equilibrium solution will then be constructed using the method of matched asymptotic expansions. It will be shown that the height of an individual spike is a function of the position of all the other spikes. The n-spike eigenvalue problem will then be derived. An n-spike solution will be shown to be meta-stable under an appropriate condition on the inhibitor diffusion coefficient. This leads to a quantization condition for the maximum number of meta-stable spikes that the system can support for a given value of £ V Finally, in Chapter 4 we will give some preliminary results for the G M system in higher spatial dimensions. In particular, we use the projection method to derive an ordinary differential equation for the location of a spike layer in a multi-dimensional setting. A variety of numerical methods and software packages were used to carry out the numeri-cal computations in this thesis. Short time simulations of the full P D E system are carried out using I M E X schemes [11, 2]. Long time simulations use the fully implicit scheme from the package P D E C O L . Numerical solutions to eigenvalue problems are computed using COLSYS and M A T L A B . 11 C h a p t e r 2 I n f i n i t e I n h i b i t o r D i f f u s i o n C o e f f i c i e n t 2.1 I n t r o d u c t i o n We now examine the Gierer Meinhardt equations in the weak coupling limit Dh —> co. We will begin by constructing a one-spike quasi-equilibrium solution. The stability of this solution will be examined by analyzing the spectrum of the eigenvalue equation resulting from a linearization about our one-spike solution. The principal eigenvalue is exponentially small and we estimate it precisely in the limit e —» 0. We then use the projection method to derive an ordinary differential equation governing the motion of the location of the spike corresponding to a one-spike solution. The case of n spikes will then be considered. The stability of an n-spike solution will be studied by a similar examination of its linearized spectrum. The scaled Gierer Meinhardt equations are given by, dP at = e2axx — a+ — , — 1 < re < 1, £ > 0 , (2.1a) am rht = Dhhxx - ph + e _ 1 — , (2.1b) ns a x ( ± l , t ) = 0, hx(±l,t) = 0. (2.1c) In the limit Dh —r oo we write h as a power series in D^1 as h = ho + D^hi + ••• . (2.2) 12 Chapter 2. Infinite Inhibitor Diffusion Coefficient Substituting this into (2.1b) we arrive at the following equations: h0xx = 0, - 1 < x < 1, (2.3a) am hixx = Thot + ph0 - C~1 — , - 1 < x < 1, (2.3b) M ± M ) = 0, (2.3c) hlx(±l,t)=0. (2.3d) From (2.3a) and (2.3c) we find that h0 = h0(t), and so h0 is spatially homogeneous. By applying a solvability condition to (2.3b) subject to (2.3d), we derive the following O D E for h0 = ho{t): 1 f1 am rh0 + ph0 - e _ 1 - / — dx = 0. (2.4) 2 y_! hs Here ho = dho/dt. We expect that the dynamics of h is much faster than that of a. Therefore, we set h0 = 0 in (2.4) and solve for the equilibrium value of ho. In this way, we get h ° = { e ~ 1 i L a m d x ) ' * 1 - ( 2 - 5 ) Thus, to leading order as Dh —> oo, the Gierer Meinhardt equations are reduced to dP at = e2axx - a+-g, - l < a ; < l , t>0, (2.6a) n0 h0= {e-l-^-j\mdxy+\ (2.6b) o x ( ± l , t ) = 0. (2.6c) This system is referred to as the Shadow System for (1.7) (see [7, 10]). To determine the range of validity of this approximation, we note that we have required hxx = 0 to be the dominant balance in a neighborhood of a spike. Thus, if we scale y = e~1(x — XQ), where x0 is the spike location, we will require that, ^ > e " 1 or Dh^e. (2.7) 13 Chapter 2. Infinite Inhibitor Diffusion Coefficient To ensure that hxx — 0 is the leading balance in the outer region, defined away from an 0(e) region near the spike, we will require that Dh 3> 1. 2.2 A O n e - S p i k e Q u a s i - E q u i l i b r i u m S o l u t i o n We now construct a one-spike quasi-equilibrium solution aE — aE(x). This solution will be symmetric about x 0 , where |x 0 | < 1, and it will achieve a global maximum at x = XQ. In addition, aE(x) —>• 0 at infinity. The quasi-equilibrium solution aE(x) satisfies e2a"E - a E + -^ = 0, (2.8a) h0 = ( e - 1 ^ J amdx^j , (2.8b) a'E(x0) = 0, (2.8c) a E ^ 0 a s i - > ±oo. (2.8d) Now we introduce the local variable y = e~1(x—x0) and we set set uc(y) = h^1 aE(xQ-\-ey), where 7 = q/(p — 1)- Substituting this into (2.8) we get the following canonical spike problem uc(y): u'c-uc + upc = 0, 0 < y < oo, (2.9a) uc —y 0 as y -» oo, (2.9b) u'c(0) = 0. (2.9c) In terms of the solution to (2.9), the quasi-equilibrium solution for (2.8) is aE(x) = hluc (e _ 1 (x - x 0 )) , (2.10a) p-i / 8\ (a + l)(p-l)-9m f°° ho=(ji) . (3 = J^umdy, J = q/(P-1). (2.10b) Here x 0 is the unknown location for the center of the spike. The existence of solutions to (2.9) can be shown by analyzing the phase plane and has been proved in [8]. 14 Chapter 2. Infinite Inhibitor Diffusion Coefficient To determine numerical values for certain asymptotic quantities below we must compute uc(y), P, and other constants numerically. To do so we first note that in the far field uc ~ ae~y as y —> oo, where a > 0 is given by (see [17]) l o g ( ^ ) , I - 1 1 log(a) P-I Jo p + i IP+I v dr). (2.11) Therefore, we can use the asymptotic boundary condition u'c + uc = 0 at y = y^, where yL is a large positive constant. To compute solutions for various values of p, we use a continuation procedure starting from the special analytical solution uc(y) — |sech 2 ( | ) , which holds when p = 2. The boundary value solver C O L N E W is then used to solve the resulting boundary value problem. In Fig. 2.1, we plot the numerically computed uc(y) when p = 2, 3,4. 1 • 6 I 1 1 1 1 1 1 1 1 r 0 2 4 6 8 10 12 14 16 18 20 y Figure 2.1: Numerical solution for uc(y) when p = 2, 3,4. We note that the solution CLE(X) will satisfy the steady-state problem corresponding to 15 Chapter 2. Infinite Inhibitor Diffusion Coefficient (2.1a), but will fail to satisfy the boundary conditions in (2.1c) by only exponentially small terms as e —> 0. This will be true for any value of x0 that is not within an 0(e) distance from the boundary. Thus, we will need to use exponentially accurate asymptotics to determine the equilibrium position of the spike. 2 . 3 T h e O n e - S p i k e L i n e a r E i g e n v a l u e P r o b l e m To examine the stability of the quasi-equilibrium spike solution found in the previous section, we will linearize about this solution and we study the spectrum of the corre-sponding eigenvalue problem. The resulting eigenvalue problem is of a non-local nature. Results from [3] suggest a numerical method for the analysis of the spectrum of such a problem. To solve the non-local problem we introduce a continuation parameter to gradually introduce the non-local effects. The eigenvalue problem on the extended real line will then be considered, for which some exact results exist. The perturbing effect of a large but finite domain will then be studied. To begin our analysis, we derive the eigenvalue problem in terms of cp and rj defined by a(x, t) = aE(x) + ext(j)(x), (2.12a) h(x, t) = h0 + extr](x). (2.12b) Here aE and h0 are given in (2.10) while <f> <C aE and t] <C h0. Substituting this into (2.1) results in the following eigenvalue problem; p—I p t2<f>xx-<f> + P^<j>-q-^iV = W, (2.13a) a m _ 1 am DhVxx ~ m + mt~1^—^ - se-1—fj-77 = rXn . (2.13b) n0 tiQ Substituting (2.10a) and (2.10b) into (2.13) we get e2(j>xx - 0 + pvP-V - qhlv-q-xvPcJ] = \<j>, (2.14a) Dhr]xx + me-lhfm-l)-su™-l(l) - s c - 1 ^ " * " 1 ^ ^ = ^Xn. (2.14b) 16 Chapter 2. Infinite Inhibitor Diffusion Coefficient We then expand r\ as a power series in D^1, r) = r}0 + D^r]l + O(Df), (2.15) and we substitute this into (2.14) and collect powers of D^1 to obtain Voxx = 0, - 1 < x < 1, (2.16a) Vixx = Wo ~ me-^^-V-'u?-1^ + se-1hT~'~1v^r}0 + rXr]0, - 1 < x < 1 (2.16b) 7 f o x ( ± l ) = 0, (2.16c) ?7ix(±l) = 0. (2.16d) Thus, ?7o is a constant independent of x. To determine rjo we apply a solvability condition on the 771 problem to get (2/x + 2s(3hlm-s-x + 2Ar)77o = e ^ m / i ^ 1 ^ f u^1^ dx , (2.17) where (3 is defined in (2.10b). Solving for 770 we get t 3 b = g m / 1 ° / (2.18) 2(/x(s + 1) + Ar) c v ' Our non-local eigenvalue problem for <f> = <f)(x), defined on |x | < 1, is obtained by substituting (2.18) into (2.14a) Lt<f> = e2cf>xx - <j> + pul-l<f> - m ^ J l q ^ l X t ) f u r 1 * * * = (2.19a) <^(±1) = 0. (2.19b) The integral term in equation (2.19) changes the nature of this problem drastically from standard Sturm-Liouville problems and non-standard techniques will be needed [3]. Also note that the eigenvalue A appears on both sides of the equation. Thus, (2.19) is an implicit eigenvalue equation. However, since r is typically very small, we may assume that r = 0 as a simplifying approximation. This approximation is used in all of the 17 Chapter 2. Infinite Inhibitor Diffusion Coefficient analysis below. However, the case of a small r would not be significantly more difficult to analyze. In (2.19) we note that uc = uc [e~l(x — x0)]- Therefore, we will only seek eigenfunctions that are localized near x — x0. These eigenfunctions are of the form 4>{y) = cf)(x0 + ey), y = e~l(x - x0) . (2.20) Therefore, we can replace the finite interval by an infinite interval in the integral in (2.19) and impose a decay condition for <j)(y) as y —> ±oo. This gives us (with r = 0) the eigenvalue problem for the infinite domain — co < y < oo Lej> =~4>yy~A K " ^ - 2Bis+l) r """^ ^ = ^' ^'^^ <j>(y) -> 0 as y -> ± o o . (2.21b) To treat the non-local eigenvalue problem, we split the operator Le into two parts, AJ> = e2ct>xx - 4> + pug" V, B<t> = Ip™^ f « r V dx. (2.22) We define a new operator Lg(j) = Acp — 5B(f>. When S — 0 we have a simple Sturm-Liouville problem. At 6 = 1 we have our full non-local eigenvalue problem. We define Lg, A and B in a similar fashion, but on the extended domain —co < y < oo with the appropriate boundary conditions at ±oo. To observe that L£ has a zero eigenvalue, we first note that if we differentiate (2.9a) with respect to y = e~l(x — xo), it is clear that Au'c = 0. In addition, uc(y) is even about y = 0 and is increasing for y < 0 and decreasing for y > 0. Thus, u'c is odd about y = 0. Therefore, J^u^^u^dy = 0, which implies that Bu'c = 0 as well. Thus, Leu'c = 0. Moreover, uc and u'c tend to zero exponentially as y —> ±oo. Therefore, the eigenvalue problem (2.21) has a zero eigenvalue with corresponding eigenfunction (j>(y) = u'c(y). Now for the finite domain problem (2.19), the function u' c[e_ 1(a; — x0)] fails to satisfy the equation and boundary conditions of this problem by exponentially small terms as 18 Chapter 2. Infinite Inhibitor Diffusion Coefficient e —>• 0. Therefore, we expect that the presence of the finite domain will perturb the zero eigenvalue and corresponding eigenfunction of the extended problem by only an exponentially small amount. The function uc(y) has a unique maximum at y = 0 and thus the eigenfunction u'c(y) has exactly one zero at y = 0. This implies that uc(y) corresponds to the second eigenfunction of A . Hence, the principal eigenvalue of A is positive and bounded away from zero. Therefore, the principal eigenvalue of A for the finite domain problem is not exponentially small. Since Lg has a positive eigenvalue when S = 0, we must consider what happens to this eigenvalue as 5 ranges from 0 to 1. If this eigenvalue remains positive then, since we expect that the eigenvalues of Lg and Lg will differ only by exponentially small amounts, we can conclude that the one-spike quasi-equilibrium solution is unstable. Alternatively, if this eigenvalue crosses through zero at some finite value of S < 1, then the principal eigenvalue of Lg when 5 = 1 (which corresponds to our eigenvalue problem (2.19) will be exponentially small. Hence, if this occurs, the one-spike solution is anticipated to be meta-stable. We now estimate an eigenvalue for the infinite domain operator Lg when 5 <C 1. To fix notation, let tpo(y) and A 0 be the first eigenpair of our local operator A and let X0(5) be the eigenvalue of Lg for which XQ(5) —> A 0 as 8 —> 0. The corresponding eigenfunction of Lg is denoted by <f>(y;5). Specifically, we will calculate the sign of A0(0) analytically. Thus, we have that tf>o(y) and <j>(y;6) satisfy <f>oyy + (-1 + P " ? - 1 ) ^ = Ao^o , (2.23a) 0o -)• 0 , as y-¥ ±oo . (2.23b) and tyy + (-1 + K _ 1 ) £ - ^BXS+I) F ^ ^ ' ( 2 - 2 4 A ) - » 0 , as y -»• ±oo . (2.24b) 19 Chapter 2. Infinite Inhibitor Diffusion Coefficient Multiply (2.23) by 4> and (2.24) by ^ 0 and subtract the resulting equations. Then, integrating this the result from — oo to oo, we arrive at the following relation = - s ^ i j £ > * j y - ^ y -Now taking the limit as 5 —>• 0, we have a% mg f ^ u ^ d y f^u^fody d6 s = 0 2P(s+l) fcfady • [ Z - Z b ) Since uc > 0 on ( — 0 0 , 0 0 ) and <f>0 is of one sign, we conclude that ^f|«j=o < 0. Thus, X0(5) — A 0 < 0 when 8 is sufficiently small. We must now examine whether this inequality, which occurs when 5 is small, will persist as 5 increases to cause A 0 to cross through zero at some value 0 < S < 1. We will now examine the eigenvalues of the non-local eigenvalue problem on the infinite line (2.21). Recall that in terms of the local and non-local operators A and B, respectively, this problem can be written as Here LS(j) = A(p- 5B(j> = Xcf) - 0 0 < ?/ < 0 0 (2.27a) 0 - 4 0 , as y-t ± 0 0 . (2.27b) = 4>w - c* + pv?-1^, B<f) = 2 ^ 1 ) f " « r ^ dy. (2.28) The calculation of the eigenvalues of this problem will require some numerical analysis. Thus, we will work with a specific parameter set. The values (p,q,m,s) = (2,1,2,0) are commonly used in simulations, so we will work exclusively with this set. For this parameter set, we begin by reviewing some exact results for the spectrum of the local eigenvalue problem A4> = X4>, (2.29a) as y - > ± o o . (2.29b) 20 Chapter 2. Infinite Inhibitor Diffusion Coefficient This problem has three isolated eigenvalues and a continuum of eigenvalues, comprising the continuous spectrum. These three isolated eigenvalues (when p = 2) are A 0 = 5/4, Ai = 0 and A 2 = —3/4 with eigenfunctions 0 O = sech2(y/2), ^ = tanh(y/2)sech2(y/2) and 4>2 = 5sech3(j//2) — 4sech(y/2), respectively (see [4]). For the corresponding finite domain problem, we note that the eigenfunctions above, written in terms of y = e~l(x — XQ), will fail to satisfy the boundary conditions in (2.19) by only exponentially small terms as e —> 0. Thus, we expect that the eigenvalues of A will be only slightly perturbed from those of A. As we have previously noted, the zero eigenvalue of (2.29) will persist as 8 ranges from zero to one. Hence, there is an eigenvalue of (2.19) that is exponentially small as e —> 0. Now we will compute the eigenvalues A0(5) and A2(<5) for which A0(5) —> 5/4 and A2(<5) —»• —3/4, as 8 —> 0. We need to compute these eigenvalues numerically. To do so, we use the initial guesses provided above for 8 = 0 and then use a continuation procedure to compute these eigenvalues as 8 increases. The computations are done using C O L N E W . In Fig. 2.2 we plot the numerically computed A0(5) and A2(<5) versus 8. As can be seen from this graph, A 0 ~ 0 for 8 = 1/2. It will be shown analytically that this relation holds with equality (i. e. A 0 = 0 when 8 = 1/2). Shortly after A 0 becomes negative, it becomes complex. At this point, C O L N E W is no longer able to track the eigenvalue. To understand what is happening, we computed A2(5) as well. Figure 2.2 illustrates the situation. As 8 increases from 0 to 1, An is decreasing and A 2 is increasing. At a value of 8 ~ 0.65 the two eigenvalues collide and split into complex conjugates eigenvalues with negative real parts. To track the eigenvalues beyond 8 w 0.65 one must employ a different numerical technique. These computations are done by discretizing the finite domain problem when e <§C 1. These eigenvalues are exponentially close to the eigenvalues of Lg and so we can neglect these exponentially small errors. The operator Lg may be approximated by a discrete linear operator (i.e. a matrix) Cg. The eigenvalues of 21 Chapter 2. Infinite Inhibitor Diffusion Coefficient the continuous problem may then be approximated by the eigenvalues of this matrix. To discretize the operator, we use the centered difference approximation of the second derivative for the local operator. The non-local operator is approximated using the Trapezoidal rule. This then results in the following matrix, o ••. 0 ^ 2 3 0 0 \ rn-2,n-3 Tn-2,n-2 Tfi-lfi-l 0 rn-l,n-2 Til-lfi-l ) ( + 8 where, r l f l = -2e2/h2, r 1 > 2 = 2e2/h2, = e2/h2, rhl = -2e2/h2 + (-1 +pupc-1((xi - x0)/e)), n,i+i = e2/h2, mquP((xi - xQ)/e) _x -uc ((-1 - x0)/e)h/2, Si,l = 2p(s + l) mqupc((xi-xQ)/e)^m_1 2p(s + l) UT ((xj ~ x0)/e)h, mqupc{(xi-x0)/e)^m_1 x. 2P(s + l) h = 2/n, - 1 + ih. urL((l-x0)/e)h/2, \ (2.30) (2.31a) (2.31b) (2.31c) (2.31d) (2.31e) (2.31f) (2.31g) (2.31h) (2.31i) (2.31J) Here n is the number of grid points. By numerically calculating the eigenvalues of £5 we give numerical results for A 0 in Table 2.1. Since the real part of A 0 remains negative as 8 —¥ 1, we conclude that the one-spike quasi-equilibrium solution is stable for this 22 Chapter 2. Infinite Inhibitor Diffusion Coefficient parameter set. Similar computations can be performed for other values of p, q, m and s. It is possible to find the critical value of 5, denoted by 8 = 8C, for which A0(5C) = 0. At this value of 5, we will have two eigenfunctions corresponding to the zero eigenvalue. One of these eigenfunctions is known to be u'c(y). Thus, we may use the method of reduction of order to find the other eigenfunction (j>(y). Introduce v(y) by (f> = vu'c. Then, in terms of v, (2.27) with A = 0 becomes u'cv" + 2v'u"c - upc5J = 0, (2.32) vu'c —> 0 as y —> ±oo . Here j _ mq ' /oo uTX<vdy. (2.33) • o o 2/3(5 + 1) We will consider I as a constant, independent of v for now. Next, we substitute w = v' in (2.32) to get the following equation for w: [W{u'cf] = 5cIvPcu'c. (2.34) The solution is 5 I up+1 C w = ^ i W + W ( 2 - 3 5 ) To satisfy the boundary conditions in (2.32) as y —> ±oo, we need only require that w is bounded as y —> ±oo. Clearly this implies that C = 0. We then have the following solution for v, SJ fy up+l s i ry P + 1 7-00 « ) pFinally, we substitute the equation above into (2.33), to obtain the following relation: 2/?(s + 23 Chapter 2. Infinite Inhibitor Diffusion Coefficient Since I ^ 0, we can cancel / to get the following expression for Sc, *•=U^t,+i)/>'* (/IS*') *)"'• (2-38) For the parameter set we used, the integral above may be evaluated exactly. Substituting Uc(y) = |sech(t//2)2, m = 2, p = 2, q = 1 and s = 0 into the equation above, results in Sc = | as is suggested by Fig. 2.2. i. 4 1 1 1 1 —1 1 1 2 -1 -0 8 0 6 0 4 0 2 0 -0 2 -0 4 -0 6 _n o 1 1 1 1 1 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 5 Figure 2.2: A 0 and A 2 versus. 5. 2 . 4 A n E x p o n e n t i a l l y S m a l l E i g e n v a l u e In the previous section, we showed that the only positive eigenvalue of the local operator A becomes negative with the inclusion of the non-local effects. Thus, for the non-local operator L e , the principal eigenvalue will be exponentially small. We denote this eigen-value by A i . To predict the dynamics of the quasi-equilibrium solution, we must obtain 24 Chapter 2. Infinite Inhibitor Diffusion Coefficient s A 0 0.0 1.2518 0.1 1.0073 0.2 0.76149 0.3 0.51345 0.4 0.26158 0.5 0.0052548 0.6 -0.28247 0.7 -.59237+ 0.15315z 0.8 -.71522+ 0.23035z 0.9 -.84093+ 0.23008i 1.0 -.98551 + 0.14507z Table 2.1: S and A 0 for the case (p, q, m, s) = (2,1,2,0). a very accurate estimate of A i . Let fa denote the eigenfunction corresponding to X\. We expect that fa ~ C\uc {e~1(x — x0)) in the outer region away from O(e) boundary layers near x = ± 1 . The behavior of fa in these regions will be analyzed using a boundary layer analysis. To begin the boundary layer analysis we write fa in the form fa(x) = C i « [e~l(x - x0)] + fa [e~l{x + 1)] + fa [ e ^ l - x)]) . (2.39) Here fa(r)) and fa(n) are boundary layer correction terms and C\ is a normalization constant given by ~*-*(/j<(y)rdy) 2 • (2.40) 25 Chapter 2. Infinite Inhibitor Diffusion Coefficient Thus, Ci = UP) , where 3=1 « ) 2 dy. (2.41) J — o o In the boundary layer region near x = —1, uc[e~l(x — xQ)] is exponentially small as e —> 0. Thus, as e —> 0, ^ ( 7 7 ) satisfies 0j' — 0£ = 0, 0 < 77 < 0 0 , (2.42a) ^(0) ~ - a e - £ _ 1 ( 1 + x o ) . (2.42b) Similarly, the boundary layer equation for 4>r(rj) is 4'r - (j>r = 0, (2.43a) <^ (0) ~ a e - e _ 1 ( 1 - X o ) . (2.43b) Here a is defined in (2.11). Solving the boundary layer equations we get ^ ( 7 7 ) = ae-e~1{1+ao>e-"> (2.44a) ^ ( 7 7 ) = - a e - 6 - ^ 1 - 1 0 ^ - " . (2.44b) To estimate Ai we first derive Lagrange's identity for (u,Lev), where (u,v) = f^uvdx. Using integration by parts we derive (v, Leu) = e 2 (uxv - vxu) |^=LX + (u, L*v), (2.45) where L*v = e2vxx - v + ur'v - 2 ^ " ^ / ' «S« <*x . (2.46) We now apply this identity to the functions u'c[e_"1(x — x0)] and </>i(x) to get « , Le<h) = -efau'XzU + (fc, L e V c ) . (2.47) 26 Chapter 2. Infinite Inhibitor Diffusion Coefficient We will now examine each of the terms in (2.47). We begin with (u'c, Lefa). The dominant contribution to this integral arises from the region near x = x0 where u'c[e~l(x — x0)] ~ ^ r . Therefore, the inner product can be estimated as (u'c,Lefa) = ^-(fa,fa), cV / -A 1/2 (2.48a) (2.48b) (2.48c) since fa is normalized. Next, to estimate —efauc\x::zl_i, we will use our asymptotic esti-mates of uc and fa. Since uc(z) ~ ae - ' 2 ' as z —> ±oo we have that u'c [ e _ 1 ( ± l — x0)] ~ ae~e 1( 1 =F X°). In addition, using the previous boundary layer results for fa we get fa(±l) ~ ^2C\ae~e 1(lI?x°). Using these results and the estimate for Ci, we get -efauXz\ ~ 2 ^ a 2 ^ e - 2 £ _ 1 ( i + * o ) + e ~ 2 e ~ 1 ( 1 ~ X o ) ^ (2.49) The only term left to examine is (fa, L*€u'c). Since u'c is a solution to the local operator, we have T*ii' — — L e U c ~ 2/?(* + l) j ^upcu'cdx, mqu] m—1 2(3(s + l) p + 1 -u\ P + i x = - l ,m—1 I)(P + I) vc Thus, the term (fa,L*u'c) is approximated by amquc 2(3(s + l)(p+l( l + x o ) _ g - ( p + l ) e ' • ' ( l - s o ) ^ _ (2.50) amq 2/3(s + l)(p + l) Ciamq 20(s + l)(p + l) Ciamq 2/?(s + l)(p + l)m -(p+^e-^l+xo) _ „ - ( p + l ) e - ( p + ^ e - ^ l + x o ) _ „ - ( p + l ) e - ( p + ^ e - 1 ( l + x o ) _ P - ( P + i(i-xo)) y1 u ^ - ^ i ^ c f a , ^ l - x o ) ^ l ) e - i ( l - x o ) ^ ^ x = l x = - l g m e 1 ( l + x o ) _ g ~ m e 1 ( 1 _ x o ) (2.51) 27 Chapter 2. Infinite Inhibitor Diffusion Coefficient Since p > 1 and m > 1 the term from equation (2.51) will be asymptotically negligible compared to the term from (2.49). Therefore, to within asymptotically negligible terms, (2.47) gives us the following asymptotic estimate for Ai as e —>• 0: Ax ~ 2a 2 /T 1 (e-^d+zo) + e -2 e - i ( i -* 0 )) . ( 2 . 5 2 ) In (2.52), a and (3 are defined in (2.11) and (2.41), respectively. This is the main result of this section. This estimate holds for p, q, m and s satisfying (1.2). As an example, we take the parameter set (p,q,m,s) = (2,1,2,0). For these values we can calculate that uc(y) = |sech 2(y/2), a = 6 and /3 = 6/5. Therefore, for a spike centered at x0 = 0 with e = 0.02 we have that A 1 « 2 ^ ( 2 e - 2 / ° 0 2 ) , ~ 0.4464091171 x 10~ 4 1. We end this section with a few remarks. Firstly, we recall that Ai and fc ~ C\uc [e~l(x — x0 are an eigenpair of Lg when 5 = 0. To within negligible exponentially small terms this eigenpair remains an eigenpair of Lg as 5 ranges from 0 to 1. To see this, we note that the only difference between the calculations of the eigenvalue for the local problem and for the non-local problem, is that the term (L*u'c, fc) in (2.47) would be replaced by (A(pi, fc) = 0, since A is self-adjoint. In the final calculation of Ai the term {L*u'c, fc) was ignored since it is asymptotically exponentially smaller than the other terms in (2.47). Secondly, we note that (Ai, fc) is an eigenpair of the adjoint operator, L*. For the same reasoning as above, fc would have the same interior behavior near x = XQ and the same boundary layer correction terms near x = ± 1 . Repeating the calculation to find A^, we would arrive at the same estimate as in (2.52). 28 Chapter 2. Infinite Inhibitor Diffusion Coefficient 2 . 5 T h e S l o w M o t i o n o f t h e S p i k e The quasi-equilibrium solution fails to satisfy the steady-state problem corresponding to (2.1) by only exponentially small terms for any value of XQ in \x$\ < 1. Moreover, the linearization about this solution admits a principal eigenvalue that is exponentially small. Therefore, we expect that the one-spike quasi-equilibrium solution evolves on an exponentially slow time-scale. We will now find an equation of motion for the center of the spike corresponding to the quasi-equilibrium solution. To do so we first linearize (2.1) about a(x,t) = / i Q U c [ e _ 1 ( x — x0(t))], where the spike location x0 = x0(t) is to be determined. For a fixed x0 we have shown that the linearization around this solution has an exponentially small principal eigenvalue as e - > 0. By eliminating the projection of the solution on the eigenfunction corresponding to this eigenvalue, we will derive an equation of motion for xo(t). This procedure is known as the projection method and has been used in other contexts (see [15], [17], [14] and [16]). To proceed with the analysis, we will need to use the orthogonality property the eigenfunc-tions. However, it is clear that the operator Le is not self adjoint, so the eigenfunctions may not be orthogonal with respect to the standard inner product. However, the local operator is self-adjoint and therefore has a complete set of orthonormal eigenfunctions. As previously noted, the principal eigenpair of Le corresponds to an eigenpair of the adjoint operator L*. Moreover, it is also the second eigenpair of the local operator A. We will refer to the eigenpairs of the local and adjoint operator as (A;,fc) and (X*,4>*), respectively. We are now ready to examine the motion of a spike. We begin by linearizing around a moving spike solution by writing, a(x, t) = CIE(X] x0(t)) + w(x, t), where CLE{X; x0(t)) = h],uc [ e - 1 ( x — x0{t)] , (2.53) 29 Chapter 2. Infinite Inhibitor Diffusion Coefficient and w <C aE. We also assume that wt <C e~lXQaE, where x0 = dx0/dt. We substitute this ansatz into (2.6) to get Lew = -e~lx0hlu'c [<Tl{x - xQ)] , -1 < x < 1, rj > 0 (2.54a) wx(-l,t) = - e - ^ X [ e _ 1 ( - l -x0(t))] , (2.54b) wx(l,t) = -e-xhlu'c [e-l{\ - x0(t)] . (2.54c) Here Le is the operator defined in (2.19). Next, we expand w(x,t) as an eigenfunction expansion in terms of the eigenfunctions of our local operator A, w(x,t) = f2^ir&- (2-55) i=i * Using the orthonormality of the eigenfunctions we may isolate the coefficient of Di, D1(t) = \1(w,<f>1), (2.56) since (A~i,fc) = (Ai,fc) . We also know that (AJ,<^*) = (Ai,fc) . Thus, we may write the expression above as D1(t) = (w,L*</>1), = - e 2 f cw x | *=^ + (L e u;,fc), = -e 2fc^xlx=l-i - e - ^ o ^ K ^ i ) . (2.57) We can calculate wx(±l,t) from (2.54). Then, using our asymptotic estimates for fc, the equation above for Di becomes £>,(*) - - e 2 d a ^ ( e - ( 1 + I o ) / V c [(1 + x0)/e] + e^ 1 "* 0 ^ [(-1 - x0)/e] - e^xoh&Cr1, ~ -2^Jla2hl (-e-2^+xo^ + e - 2 ( i - x o ) / e ) _ -y/|^ X o . (2.58) Since A _ 1 L>ifc is 0(1) in the region near x = x0, we must impose the solvability condition that D\(t) —)• 0 as e —> 0. Setting £>! = 0 in (2.58) yields the following equation of motion 30 Chapter 2. Infinite Inhibitor Diffusion Coefficient for the center of the spike x0 = x0(t): x0(t) ~ ^j- [e-^+x°y< - e - 2 ( i - o ) A ] . ( 2 . 5 9 ) This is the main result of this section. Setting x0 = 0 we find the equilibrium position of the spike to be located at x0 = 0 and it is stable. 2 . 6 A n n - S p i k e S o l u t i o n We will now examine the properties of an n-spike quasi-equilibrium solution. The anal-ysis will proceed in the same manner as for the case of the one-spike quasi-equilibrium solution. The stability of an n-spike quasi-equilibrium solution will be examined by linearizing about this solution and studying the resulting spectrum. We begin by defining an n-spike quasi-equilibrium solution by n - l anMx) = hl,E^2uc - xt)] , (2.60a) t = 0 K,E = (e- 1 ^ J1 alE dx^j ^ , (2.60b) where 7 — q/(p — 1). Substituting (2.60a) into (2.60b), we can determine hn>E as ( n0\ ( a + l ) £ - l ) - , m M ) ' ( 2 ' 6 1 ) where @ was defined in (2.10b). In (2.60a), the spike locations Xi for i = 0,.., n — 1 satisfy — 1 < x0 < xi,.., < xn-i < 1. They correspond to local maxima of a n ,£ . We now linearize (2.1) about an>E and hn>E by introducing </> and n defined by a(x, t) = antE(x) + ext(j>(x), (2.62a) h(x,t) =hn,E + extn(x). (2.62b) (2.62c) 31 Chapter 2. Infinite Inhibitor Diffusion Coefficient Here (j) ^ an,E a n d rj <C hn>E. Substituting (2.62) into (2.1) we get the following eigen-value equation e2fax ~ & + p-jj^(t> ~q-r^V = AcA, (2.63a) nn,E h n E am~l am DhVxx ~ m + m e _ 1 T J ^ - ^ - se~1-^[rj = Arn. (2.63b) nn,E hn,E Since each spike of the quasi-equilibrium solution is localized to within an 0(e) region near x — Xi for some i, we look for an eigenfunction fax) of the form 71-1 fax) = YJ&V~\x - ^)] • (2.64) i=0 Therefore, we need to introduce local coordinates near each spike. In particular, the ith set of inner variables are defined as fa(Vi) = Hxi + eVt), Vi = e~l(x ~ xi)- (2-65) Once again, we expand 77 as a power series in D^1, V = m + D-1m+0(D-2). (2.66) Substituting this expansion into equation (2.63) we get the following equations for 770 and Voxx = 0 , - K x < l (2.67a) (CZ1 , a™ mxx = prjo - me '-^—(1) +se-'—^rjo + rXrjo, - 1 < x < 1, (2.67b) nn,E h n E Vox(±l) = 0, (2.67c) mx(±l) = 0. (2.67d) Thus, 770 is a constant and it can be determined by imposing a solvability condition on the problem for 771. This condition requires that /: ( /irjo-me x-^-</> + se 1T^VO + rXr]0 ) dx = 0 . (2.68) 32 Chapter 2. Infinite Inhibitor Diffusion Coefficient The integral is decomposed into the sum of four separate integrals. We then can calculate the third integral as -5^odx = K m E s - % / Wn)dVi 1 nn,E j = 0 • ' - o o = 2nrj0sPhlmE-s-1. (2.69) Substituting (2.69) and (2.64) into (2.68) we can determine rjo as mh^m~^~s n~l r°° % = — T T T ^ ^ I E ur1(yl)Uyi)dyl. (2.70) Substituting (2.70) into (2.63a), we arrive, after a lengthy algebraic calculation, at the following eigenvalue problem corresponding to an n-spike solution: p 71—1 poo <?<\>xx + pupc4> ~ 2n/3(/x(s + 1) + rA) ^ J U ^ ~ L D ^ = A<^> \ X \ < 1 (2.71a) fc(±l)=0. (2.71b) Finally, we use localized coordinates to examine the stability of each spike. This yields on the interval — oo < yi < oo that 71—1 p Tl — 1 /»0O fam ~4> + PY1UC le~l (x ~ $ ~ 2n^( / i (sT 1) + T A) ^ / uTl{Vi)^i{yi)dVi = A^> M < c (2.72a) $ -»• 0 as ?/, ->• ±oo . (2.72b) Since r is typically very small, we can set r = 0 in (2.72) as a simplifying approximation. Now we note that if each fc were independent of i (i. e. fc(?/i) = $(2/i)) f ° r « = 0,.., n — 1, then S^To J"^ u ™ - 1 ^ ) ^ ) % = n u™" 1( ?/)$(y) dy. The factor of n would cancel in (2.72) and we would be left with the same eigenvalue problem as (2.19). Thus, for the parameter set we have used previously, we would conclude that an n-spike solution is 33 Chapter 2. Infinite Inhibitor Diffusion Coefficient meta-stable. However, we now show that this conclusion of meta-stability is erroneous. To see this we note that we can construct a global eigenfunction by taking 4>i(yi) = b^(y) for some constant b{. The non-local term in (2.72) then becomes "—1 poo poo I 1 \ £ / <~\yi)Uvi)dyi = / u™-\yMy)dy J > . (2.73) i=0 J~°° J-°° \i=0 J Then, if we impose the constraint that 71-1 X > = 0 , ( 2 - 7 4 ) 7=0 the non-local term vanishes. Hence, with this constraint, $(?/) satisfies the local eigen-value problem $" - $ + pup-l$ = A 0 $ . (2.75) This problem has exactly one positive eigenvalue A 0 . When p = 2, we found that A 0 = 5/4 with corresponding eigenfunction $o(y) = sech2(y/2). Hence, under the constraint (2.74), A 0 is also a positive eigenvalue of (2.72). This then leads to an instability. In summary, when there is more than one spike we may always construct an eigenfunction of the form <j>(x) = YJlZi b& [e~\ x — Xi)] where YLi=o ^ = 0- This eigenfunction has a positive eigenvalue. Therefore, it is impossible to find a stable multiple spike solution for large values of Dh. We now illustrate this instability result numerically for a two-spike solution for the pa-rameter set (p,q,m,s) = (2,1,2,0), fx = 1, r = 0.01, Dh — 40, and e = 0.05. We took the quasi-equilibrium solution as our initial condition. The first spike (Spike 1) is centered at x0 = —0.5 while the second spike (Spike 2) is centered at x\ = 0.5. In Table 2.2 we tabulate the numerically computed amplitudes of the two spikes as a function of time. We now use this data to estimate the positive eigenvalue. We remark that the data in Table 2.2 is taken after the simulation has been run approximately t = 20 units 34 Chapter 2. Infinite Inhibitor Diffusion Coefficient to eliminate any transients and to ensure that the positive eigenvalue is dominant. After this time the solution at the spike locations x — x0 and x = X\ will be approximately given by, a{xl,t) ^ a2,E{xi) + eXotfa{xi), i = 0 ,1. (2.76) This relation will only govern the linear instability of CL2,E- For the parameter set we have used a2tE{xi) = 6.25. Then, we can re-write (2.76) as, A 0 t + log[^o(a; i ) ]«log( |a(a; i ,*)-6.25 | ) , i = 0 ,1. (2.77) To estimate A 0 from the data in Table 2.2, we take x\ — 0.5 and evaluate (2.77) at two different values of time, labeled by t\ and t2. Using the numerically computed values for a(0.5,t) at t = ti and t = t2 gives us two equations for the two unknowns 0o(O.5) and A 0 . In this way, A 0 can be estimated. In Table 2.3 we give the numerical results for A 0 and 0o(0.5) using various values of t\ and t2. For this parameter set, we would expect that the principal eigenvalue is 1.25. The interpolated values, obtained by our numerical procedure, are all close to 1.25 as expected. 35 Chapter 2. ' Infinite Inhibitor Diffusion Coefficient time Spike 1 Height Spike 2 Height 19.5 6.2738663390032 6.2635545772640 19.8 6.2761841264723 6.2612374542097 20.1 6.2795439171902 6.2578790593718 20.4 6.2844142872978 6.2530116219365 20.7 6.2914746492378 6.2459574213615 21.0 6.3017102226534 6.2357347923334 21.3 6.3165498360813 6.2209223724487 21.6 6.3380658088900 6.1994635220925 21.9 6.3692634678024 6.1683858288480 22.2 6.4144988374999 6.1234022942033 22.5 6.4800753144382 6.0583539884737 22.8 6.5750761790975 5.9644587136514 23.1 6.7124619645111 5.8293778760449 23.4 6.9103141671528 5.6362910167926 23.7 7.1926098041313 5.3636802602846 24.0 7.5876099073842 4.9877041819062 24.3 8.1196649411629 4.4906888214834 24.6 8.7900467006609 3.8780493118962 24.9 9.5540932480348 3.1943794842134 25.2 10.3232394145038 2.5150430348060 25.5 11.006159488840 1.9098035904776 Table 2.2: Height of spike 1 centered at x0 = —0.5 and of spike 2 centered at x\ — 0.5. 36 Chapter 2. Infinite Inhibitor Diffusion Coefficient h * 2 a(.5,*i) a(.5,t2) Ao 0o(-5) 22.8 23.4 5.9644587136514 5.6362910167926 1.275223721 _ e - 3 0 . 3 2 8 4 6 9 4 9 23.1 23.7 5.8293778760449 5.3636802602846 1.242238171 _ e - 2 9 . 5 6 1 7 2 2 1 5 22.5 23.7 6.0583539884737 5.3636802602846 1.276189822 _ e - 3 0 . 3 6 6 3 7 6 2 9 22.2 23.4 6.1234022942033 5.6362910167926 1.315422050 _ e - 3 1 . 2 6 9 1 1 0 3 8 Table 2.3: Logarithmic Interpolation of A 0 and </>n(.5). 37 C h a p t e r 3 Finite Inhibitor Diffusion Coefficient In the previous chapter, we examined the Gierer Mienhardt equations in the limit e —> 0 and Dh —> oo. In this chapter, we analyze the case of a finite Dh in the limit e —>• 0. From previous numerical experiments, it would seem that a smaller inhibitor diffusion coefficient can lead to more spikes that are stable. We begin with the scaled Gierer Mienhardt system (see (2.1)) from the previous chapter, ap at = e2axx — a + —, — 1 < x < 1 =,, rj > 0, (3.1a) hq am rht = Dhhxx- nh + e _ 1 —, (3.1b) hb where p,q,n,s satisfy (1.2). We construct a quasi-equilibrium solution to (3.1) with n spikes using the method of matched aysmptotics. To examine the stability of this solution, we study the associated eigenvalue problem arising from linearizing (3.1) about our quasi-equilibrium solution. A n inner solution in an 0(e) neighbourhood of each spike is matched to an outer solution defined away from the spike. The cases of one spike and of n spikes(n > 1) will be treated separately. 3.1 A O n e - S p i k e Q u a s i - E q u i l i b r i u m S o l u t i o n In the limit e —> 0, we construct a quasi-equilibrium solution to (3.1) with exactly one spike. The spike is centered at x 0 , with — 1 < x0 < 1 and x0 is taken to be the local 38 Chapter 3. Finite Inhibitor Diffusion Coefficient maximum of a. We use the method of matched asymptotics to construct the quasi-equilibrium solution. In the inner region, defined in an 0(e) neighbourhood of x0, we introduce the following inner variables, a(y) = a(x0 + ey), h(y) = h(x0 + ey), y = e _ 1 (x - x 0 ) . (3.2) Substituting (3.2) into (3.1) results in the following inner equations, aP ayy — a + — = 0, —oo < y < oo, (3.3a) hfl am Dhhyy — e2ph + e— = 0, —oo < y < oo. (3.3b) hs We then expand h and a in powers of e, h = ho + eh1-\ , a = a 0 + O(e). (3.4) Substituting (3.4) into (3.3) and collecting powers of e, we find, h0yy = 0, —oo < y < oo, (3.5a) 1 am hyy = -oo < y < oo. (3.5b) To match to the outer solution constructed below, we will require that h0 does not grow linearly in y as y —>• ±oo. Thus h0 is a constant independent of y. Therefore, So satisfies (3.3a) with an unknown, but constant value of h, i . e. h ~ h0, Thus, as in (2.10a) the quasi-equilibrium solution aE{x) = do is aE{x) = hluc [e _ 1(x - x0)] 7 = q/{p - 1). (3.6) Here uc(y) is the canonical spike solution satisfying (2.9a). To determine h0 we must match the inner solution to the outer solution, which we will construct below. To obtain a matching condition for the outer solution, we integrate 39 Chapter 3. Finite Inhibitor Diffusion Coefficient (3.5b) from y = — oo to y = oo to get 1 . r°° h[(oo) - h[(-oo) = -—hlm-s / < ( y ) dy. (3.7) Dh 7-00 Since ho is a constant, we get from (3.4) that, h'(oo) - h'(-oo) = -—hl m ~ s / um(y) dy + 0(e 2). (3.8) Dh J-OO Now, we construct the outer solution defined away from an 0(e) neighbourhood of x = xo-Since a is exponentially localized to an 0(e) region about x0, we get to within negligible exponentially small terms, that a = 0 in the outer region. In the outer region we get, to within exponentially small terms, h satisfies Dhhxx — ph = 0 on [—1,1] subject to continuity and jump conditions that must hold at x = XQ. TO derive these conditions we write the matching condition between the inner and outer solution as, h(x) ~ ho + ehi(y) + • • • , as x —> XQ, y —> oo, (3.9a) h(x) ~ h0 + dii(y) -\ , as x —> XQ , y —> —oo. (3.9b) Therefore h(x0) ~ h0. Now by subtracting (3.9a) from (3.9b) and substituting in (3.7), we then get the jump condition, 1 roo 1 7m—s /o o u?(y)dy. (3.10) •00 Where [v] = V(XQ) — V(XQ). Therefore, the outer approximation for h satisfies, Dhhxx — u.h = 0, —l<x<l, (3.11a) hx{-l) = Ml) = 0> (3-Hb) [h] = 0, (3.11c) W] = -^rhlm-s ^ u™(y)dy. (3.11d) l>h J-00 The solution to (3.11) is given by, * = — = 2 ^ = c o s h ( J^-(a;< + l ) ) c o s h ( , / - ^ - 1)). (3.12) V ^ s i n h ( 2 , / f ) V A V A i 40 Chapter 3. Finite Inhibitor Diffusion Coefficient Where (5 is defined in (2.10b), £< = mm(xo,x) and x> = max(x 0, x). In terms of this solution, ho is given by h0 = h(xo). Thus, using (3.12) we get, 1+s yn ho = I ^ c o s h ( , / £ ( x 0 + l ) ) c o s h ( A / ^ ( x 0 - 1)) I . (3.13) VL^s inh (20y y D h * D h In summary, the one-spike quasi-equilibrium solution is given by, aE = hluc(e 1(x - x0)), (3.14a) hE = 12 cosh( J-fr{x< + l))coshU-^-(x> - 1)). (3.14b) See figures 3.1, 3.2 and 3.3 for plots of outer solutions. 0.4 1 1 1 1 1 1 1 1 1 r 0.15 h 0.1 \-0.05 - | i i 0 I I I I I ^1 i lis I I I I -1 -0 .8 -0 .6 -0 .4 -0 .2 0 0.2 0.4 0.6 0.8 1 Figure 3.1: Outer Solutions for Dh = 1. We close this section with a few remarks. Firstly, we can re-establish our approximate formula for h0 when Dh ~> 1 using our formula for h0. We examine the limit as Dh tends 41 Chapter 3. Finite Inhibitor Diffusion Coefficient 0.14 0.1 -1 -0 .8 -0 .6 -0 .4 -0 .2 0 0.2 0.4 0.6 0 Figure 3.2: Outer Solutions for Dh = .1. to oo and we find, lim j — — — — — -==-D^°° \VDhl2sinh(2J^-) c o s h ( V D~h(x°+ 1))cosh(y^;(a:o _ 1+5—yrn 28\ (3 + l ) ( p - l ) - ? T > 2^ J (3.15) Which corresponds with the value found in the previous chapter(the length of our interval is 2). Figures 3.4 and 3.5 illustrate the behaviour of ho as Dh and x0 are varied. Secondly, the derivation leading to the jump condition (3.10) can be significantly short-ened by making the following observation. In the outer variables, a is localized to within an 0(e) neighbourhood near x = x0. In the inner region, h ~ h0, where ho is a constant. Therefore, for the outer equation for h, the term e - 1 a m / / r s in (3.1) has the effect of a multiple of a delta function centered at XQ as e —>• 0. To find the multiple of 5(x — x0) 42 Chapter 3. Finite Inhibitor Diffusion Coefficient 0.5 0.35 0.3 0 .25 0.2 0.15 0.1 0.05 -1 I I I T" h 1,E _J I L. ~i 1 1 r a\,E _i i_ -1 -0 .8 -0 .6 -0 .4 -0 .2 0 0.2 0.4 0.6 0. Figure 3.3: Outer Solutions for Dh = 10. (3.16) we integrate e 1am/hs over a small neighbourhood centered at x0, lim / —— dx = h1™'8 / uc(y) dy = hln~s2(3. s^°Jx0-s e/i s 7 - o o Therefore the term e _ 1 a m /h s in (3.1) may be replaced by h0rm~s2P5(x—x0) in constructing the outer solution. This yields that the outer approximation to h satisfies, hxx - ph = -h1™ s2(3S(x - x0), - 1 < x < 1 M - i ) = hx(l) = 0. It is clear that system (3.17) and system (3.11) are equivalent. (3.17a) (3.17b) In the derivation above, we have required hxx — 0 to be the dominant balance in the neighborhood of a spike. We will thus require that condition (2.7) still hold. Away from a spike, we have no longer assumed that hxx = 0 and thus the condition that Dh 3> 1 is 43 Chapter 3. Finite Inhibitor Diffusion Coefficient Figure 3.4: h0 versus Dh, for x0 = 0, ± .5 . no longer required. It is important to emphasize that (3.14) satisfies (3.1) up to exponentially small terms for any x0 E (—1,1). It also fails to satisfy the no-flux boundary conditions at x = ± 1 by only exponentially small terms for any x0 G (—1,1). Determining x0 requires expo-nential precision in the asymptotic solution to eliminate this near translation invariance. This also suggests that the quasi-equilibrium solution could be meta-stable. In order to examine these issues, we will examine the spectrum of (3.1) linearized about (3.14). 44 Chapter 3. Finite Inhibitor Diffusion Coefficient 0.3 0.25 ha °-2 0.15 1 1 1 1 — i — i — i — i — Dh = o.i Dh = l -Dh = 10 / i i i i -0.8 -0 .6 -0.4 -0.2 0 X0 0.2 0.4 0.6 0. Figure 3.5: h0 versus XQ for Dh = .01,1,10. 3.2 A O n e - S p i k e E i g e n v a l u e P r o b l e m To examine the stability and the dynamic properties of the quasi-equilibrium solution constructed in the previous section, we now analyze the spectrum of the operator derived by linearizing (3.1) about our quasi-equilibrium solution. We thus define, a(x, t) = aE(x) + ext(j)(x), h(x,t) = hE(x) + extr)(x), (3.18a) (3.18b) 45 Chapter 3. Finite Inhibitor Diffusion Coefficient where ^ « and n -C hE. We substitute (3.18) in (3.1) and linearize to get the following eigenvalue problem: p—I p e2fax - 0 + ^ f - 0 - q - ^ V = A0, (3.19a) a m _ 1 a m D h r j x x - w + n-^—<j)-S6-1-£Ir} = Xrn. (3.19b) enE n E Again, we will use the method of matched asymptotics to match an inner region about the spike to an outer region away from the spike. A similar matching condition to (3.10) will be derived and used to solve for n in terms of 0. This then leads to a non-local eigenvalue problem analogous to (2.19). We define the following inner variables by fj(y) = n(x0 + ey), fay) = (j)(x0 + ey), y = e ~ 1 ( x - x 0 ) . (3.20) Substituting (3.20) in (3.1) and noting that, in a neighbourhood of XQ, hE ~ ho and a E ~ tiluc{y), we arrive at our inner eigenvalue problem, fay - 0 + pujr 1 ^ - <M-lupcfi = Xfa (3.21a) ^Vyy - M + j W ^ - ' r C - 1 * ~ i K T - ' - ^ c V = Arr}. (3.21b) Expand fj and 0 in an e power series, fj = fj0 + efj2 + . . . , 0 = 0o + 0(e). (3.22) Substituting this into (3.21) give us the following 0(e~2) and 0(e~1) equations on —oo < y < oo: fjoyy = 0, (3.23a) Dhf)lyy + m h l { m - l ) - s u ™ - l f a - s h r - ^ f i o = 0. (3.23b) 46 Chapter 3. Finite Inhibitor Diffusion Coefficient As was the case for constructing the quasi-equilibrium solution, in order to match to the outer solution found below, we will need to eliminate the linear growth in 770. Thus, we have that 770 is a constant independent of y, which will be determined by matching to the outer solution. By integrating the 0(e~1) equation we get the following condition, which will be used just as in (3.9) to provide a jump condition for the outer equation: 77^(00) - T M - O O ) = ~ (^hT-'^fh /_°° < dy - mhtm~l)-s £ v£-% dy) . (3.24) The situation here parallels exactly the analysis leading to (3.11). By following the same matching procedure the equation above will lead to the following jump conditions for 77 at x — xn : [77] = 0, (3.25a) [Vx] = lJh (shom~S~1^ J°° <dy - mht1^-* f_ um~% dy) • (3.25b) In the outer region, we have that a# is exponentially small. Therefore, by applying the jump conditions written above, we get the outer problem for 77 Dhr]xx - m = T A«, - 1 < x < 1, (3.26a) % ( - l ) = Tfe(l) = 0, (3.26b) [77] = 0, (3.26c) [Vx] = (fhr~'-% y00 < dy - mhtm-l)-s £ < - ^ 0 dy) . (3.26d) In the outer region (j) = o(er) for any r > 0.. with <f> = 0, at all orders of e. Solving (3.26a) results in the following, 77 = A cosh cosh O^p-) , (3-27) 47 Vo -l + (shlm-s-L2pJ-oo where a c o s h ( s ^ ) c o s h ( s £ ± ) A,sinh(J) a Chapter 3. Finite Inhibitor Diffusion Coefficient where A = n (mhom~l)~S dV ~ shT'^Vo f <dy _ _ D h sinh(-) V 7 - o o 7 - o o / ]] U- + T\ (3.28) We can determine r}o = "(^o) by evaluating (3.27) at x = XQ to get , » . 7 ( n * - l ) - s / • o o - i ; —4 /"OO 3 ^ , ^ — / ur^dy, (3-29) 2D J - O O C = V , ' ^ * : A • (3-30) Finally, we arrive at our eigenvalue problem by replacing 77 in (3.21a) by (3.29) to get the inner eigenvalue problem <f>0yy + (-1 + ur'Uo - < , \ ° l n a / « r V o d y = Afc , |yI < 0 0 , (3.31a) 1 + Csn0' 2p J-00 0->O as y - > - ± o o . (3.31b) In (3.31) we note that £ depends on a, which in turn depends on A Hence, the eigenvalue is implicitly defined by (3.31). However, since r is taken to be small, we will take r = 0 to get a = (Dh/p)1/2. Ignoring the implicit A terms will result in a small error. For a more accurate results with non-zero r, an iterative approach may be used starting with the guess found by assuming r is zero. If we compare (3.31) with the eigenvalue problem (3.31) found in the previous chapter, we note that the two problems are very similar. The only difference is the coefficient that multiplies the non-local part of the operator. The quasi-equilibrium solution will thus be stable if the coefficient of the non-local component of the operator is large enough to cause the first eigenvalue to become negative. In the previous chapter we examined the 48 Chapter 3. Finite Inhibitor Diffusion Coefficient eigenvalues of our system for the parameter set (p,q,m,s) = (2,1,2,0) and found that the eigenvalue problem u2 f°° L(f> = <!>„, + (-1 + 2uc)<j> -8-± uc<f>dy = \<f>, (3.32) J-oo with cj) —>• 0 as y —>• ±oo has a zero principal eigenvalue when 8 > \. With this parameter set, (3.31) becomes 2CJIQU2 f°° (j)oyy + (-1 + 2uc)0o T^jr1 / uc<j)0 dy = A0o- (3.33) Thus, the first eigenvalue of this problem will be zero when 2(h0 > | . If we substitute the values of values of ( and h0, given in (3.30) and (3.13), into this inequality, we conclude that one spike will be stable when 2(5 > \. For this parameter set, ft = 3. Therefore, the one-spike solution will be stable for all values of Dh- As Dh —> oo this result agrees with with the results of the previous chapter. We note that the eigenvalue problem found for the case of infinite inhibitor diffusion may be re-derived by examining (3.31) in the limit as Dh —> oo. As Dh - ) o o w e have, lim C = —^-—. (3.34) Dh^oo^ 2(/i + rA) Substituting (3.34) and (3.15) into (3.31) simplifies to the result found in the previous chapter (see (2.19)). 3.3 A n n - S p i k e S o l u t i o n We now construct an n-spike quasi-equilibrium solution to (3.1). The spikes are centered at Xi for i = 0 . . . n — 1, where x0 > —1, xi+i > Xi and x n _ i < 1. We begin by defining the following sets of inner variables near each x^. Oi(yi) = a(xi + eyt), hi(y) = h(xi + eyt), y{ = e'1 (x - x{). (3.35) 49 Chapter 3. Finite Inhibitor Diffusion Coefficient Our ith inner equations, defined on \yi\ < oo, are now aiyy - dj + -± = 0, (3.36a) K Dh - 1 a™ —hiyy-phi + --^ = 0. (3.36b) We expand hi and d, as a power series in e, ^ = ^ 0 + 6 ^ 1 + . . . , di = d i 0 + O(e). (3.37) Collecting powers of e produces the following equations on < oo /lioro = 0, (3.38a) 1 d m hilyy = - — ^ . (3.38b) To match to the outer solution, as before, we will need to eliminate the linear growth in hio as y —> ±oo. We thus have that hio is constant independent of y and clearly dio = h]0uc with uc as defined previously. It is possible to use matching of the inner and outer regions to find the jump conditions, which in turn will result in a system of equations for the h^s. However, it is less awkward to use the derivation from (3.17). In the outer region will behave like 2Ph]Qm~s5(x — x{) about each Xj, as was shown in (3.17). Matching the inner and outer solutions of the h equation will thus be equivalent to solving the following problem n - l Dhhxx -ph + 2pY^hla~SKx ~ xi) = 0, - 1 < x < 1, (3.39a) i = 0 hx(-l) = hx(l) = 0, (3.39b) h{xt+) = h(xi-), (3.39c) where h^ = h(xi) for % = 0,.., n — 1 are to be determined. Since using these equations are notationally simpler we will use (3.39) to determine the values of hio for i = 0,.., n — 1. 50 Chapter 3. Finite Inhibitor Diffusion Coefficient The solution to (3.39) on the interval (x;, xi+i) for 0 < i < n — 2 is ^ sinh (y (^xi+i - x)) ^ sinh (^ /^(a ; - x^)) / / \ i ' " l - M , U / , x • s i n h v v ik (Xi+i ~Xi^) s i n h v v D~SXi+i ~x^) In the intervals near the endpoints, h is given by cosh(y^(*+i)) (3.40) h = hQ, cosh { ^ ( x 0 + 1)) ' -1 < x < x 0 , c o s h ( ^ ( l - x ) ) h = hnfi / — , x„_i < x < 1. (3.41a) (3.41b) cosh (A/7^! ~ Z n - l ) ) By integrating (3.39a) across each Xj, we get the following jump condition at each X;, hx(xi+) - hx(Xi-) =-^-h]Qm-s• (3.42) Applying (3.42) at each x;, i = 0 . . . n — 1 yields the following nonlinear system for the h •i,0-( ocu &12 0 a2i a22 a 2 3 o ••. 0 V : ' • « n - 2 , n - 3 « n - 2 , n - 2 « n - 2 , n - l 0 0 & n _ i , n _ 2 Q ! n _ i ) n _ i J , (3.43) 51 Chapter 3. Finite Inhibitor Diffusion Coefficient where the coefficients in this matrix are defined by ftn = coth( ) + tanhf J, a a « 1 , 2 M X\ ~ %0 X a —csch( ), a ctij = coth( ) + coth( ), a a c i n - l , n - 2 &n-\,n-i = tanh( ), a M Xn—1 ^ n - 2 \ a ) + C0th( ), a (3.44a) (3.44b) (3.44c) (3.44d) (3.44e) (3.44f) (3.44g) (3.44h) In general we have to solve this system numerically. However, once we have a solution to this system we can define our n-spike quasi-equilibrium solution as n - l aN,E(X) = ^2hlouc [e l(x ~ xo)] , i = i (3.45a) cosh ( 0 , 0 cosh^y^(x0+l)) ' — 1 < X < X0, K,E{.X) = { sinh( J~^{xi+i-x)\ „ sinh( ^ ^ (x - X i ) ) hio 7 ^ — <- + / i i + i , o 7~=z ^ y , Xi<x<xi+l, i = 0 . . . n -S i n h( V D^(Xi+1~XiH sinhl ^ ^-(xi+i-xOJ cosh Xn-\ < X < 1. cosh^v/^(l-xn_1)J (3.45b) We emphasize that the height of each spike will be different and will depend on the location of all of the spikes in an intricate manner. A sample calculation for 3 equally spaced spikes is illustrated in figure 3.6. 52 Chapter 3. Finite Inhibitor Diffusion Coefficient 0.18 0.14 0.12 0.1 0.08 h 0.06 Figure 3.6: A three-spike outer solution when Dh = 1. The case of infinite Dh. may again be re-derived by examining (3.39) as Dh tends to oo. To proceed with this we write h as a power series expansion in D^1, h = h0 + —hi -\ . Dh Substituting this into (3.39) results in the following equations, hoxx — 0, hixx - ph0 + 2pJ2 f^-'Six - x^ = 0, n - l i = 0 M±i) = o» M±i) = o, h0(xi+) = h0(xi-) = hiQ (3.46) (3.47a) (3.47b) (3.47c) (3.47d) (3.47e) 53 Chapter 3. Finite Inhibitor Diffusion Coefficient Thus, we have that ho is constant and hi0 = ho for i = 0 . . . n — 1. We may now use a solvability condition on (3.47b) to find h0, /l n - l (pLho -2pJ2 hlm~s5{x - Xi)) dx = 0. (3.48) This gives, n2(3\ fc°=UfJ • (3-49> which agrees with the result from the previous chapter. To find numerical solutions to the nonlinear system (3.43), we start with a large value of Dh with the initial guess from (3.49) and then use a continuation procedure on Dh to the desired value. 3.4 The n Spike Linearized Eigenvalue Problem In the limit of large Dh, it was found that it was impossible to have a stable solution with more then one spike. With Dh = 0(1), numerical evidence leads us to believe that it should be possible to find stable multi-spiked solutions. The spectrum of (3.1) linearized about an n-spike solution should confirm this. Most of the previous analysis of the single spike linearization will not change for the case of n-spikes. First we linearize about an n-spike solution by writing, a(x, t) = antE(x) + extfax), (3.50a) h(x, t) = hn,E(x) + extr)(x). (3.50b) Since anyE is exponentially small outside of an 0(e) neighbourhood of each Xi we assume that 0 may be written as 4>(x) = Y^=o 4>i(x) where each fa is localized in an 0(e) neighbourhood of Xi. In a typical inner region, again we have an,E ~ hj0uc and h n , E ~ hi0-We thus define the ith set of inner variables to be, fa(Vi) = fa(xi + eyi), fji(yi) = n(xi + eyi), = e~l(x - x{). (3.51) 54 Chapter 3. Finite Inhibitor Diffusion Coefficient Our ith inner equation in a neighbourhood of X; is thus, 4>iyy ~4>i+ pup~l^i - qhJQ~lupcfn = Afc, (3.52a) ^f)iyy - prjt + ^ ^ - ^ - u ™ - 1 ^ - -^hir^TVi = A T * . (3.52b) Again, our goal is to match the inner and outer solutions and express each pair of coupled eigenvalue equations as a single non-local eigenvalue equation for fc. We expand fji and fc as a power series in e, Vi = Vio + + • • • , fc = fco + 0(e). (3.53) Collecting powers of e we arrive at the following equations: ViOyy = 0 , (3.54a) Dhr)ilyy + mhjt~1]~S<-lko ~ sh}™-S~l<Vio = 0 • (3.54b) Again 77*0 must be a constant to match with the outer solution, which will be determined later. As in the case of one spike, we will integrate (3.54b) to obtain, W o o ) - % „ ( - ( » ) = ^ - [shr-s-lVio2(3 - mhjt~l)~S £ u^ko dy) . (3.55) Following the analysis of (3.9) we can use (3.55) to provide a jump condition for the outer solution at X,. Thus, the outer problem for 77 derived away from an 0(e) region near each spike, is DhVxx ~m = rXv (3.5 6 a ) n—l , + E UhlT'-'Wfto - mh]t-X)-s r ur%o dy) 5(x - Xi), i=0 ^ J-oo ) Vx(l) = Vx(-l) = 0. ( 3 5 6 b ) As before, we can solve for fji on each sub-interval and then apply the jump condition at each X; to arrive at a linear system of equations for {fjio}™=i- Proceeding with this 55 Chapter 3. Finite Inhibitor Diffusion Coefficient analysis results in the following tridiagonal system; «11 «12 CX.21 C*22 o ••. 0 «23 0 0 ^ I ».,0 ^ v where, ^71-2,71-3 C*„_2,n-2 &n-2,n-l 0 o; n_i ) n_2 a n _ i i n _ i J \ Vn-1,0 J \ J = - coth( ) + tanh(— ) -\ 5 tL • 1 .xx - x0 «i,2 = csch( ), a a a u - i = csch( ——^ — h , ay = I (coth(5±i^i) + c o t h ( £ i ^ = i ) N ) + ^ r - ' W c*i,i+i — —csch( ), a a On-l.n-2 = CSChf ), a a C^n—l,n—l — - ( t a n h ( ^ - ^ - ) + c o t h ( ^ i ^ ^ ) > ) + s K - ~ ^ W , mhl 7 ( m - l ) - s » o o •iO Dh /oo oo fao dy. (3.57) (3.58a) (3.58b) (3.58c) (3.58d) (3.58e) (3.58f) (3.58g) (3.58h) Let G = A'1 where A — (a^) is the matrix in (3.57). Note that when the spikes are equally spaced , this matrix is strictly diagonally dominant and is thus invertible. In terms of the entries Ojj of ©, we can solve for the r^o's in (3.57) to get, 71-1 mh 7(771—1)—s •30 Dh u™ l(t>j0dy i = 0 . . . n - 1 . (3.59) Here © y will depend on all of the Xj ' s , on Dh, on \i and on r . About each spike, our local eigenvalue problem is given by, 71-1 lh J. /»00 faoyy ~ fao + pup~lfaa - -r-hj^ul V hfQm~1]~s / u™-1^ dy = Xfa0. uh r-f 7-oo (3.60) 56 Chapter 3. Finite Inhibitor Diffusion Coefficient Each local eigenvalue equation is coupled to all other local eigenvalue problems. To use the analysis of the previous section, we need to uncouple this system. To accomplish this, we express (3.60) in matrix notation. We thus write, 00 \ (h h = \ 4>n-l J Equation (3.60) can now be written as, 0 \ (3.61) hn-l J 4>yy - 4> + Pul~14> - ^ f r - 1 e h * m - 1 > - X r <-^dy = A / 0 , Uh J —oo (3.62) where 0 = (Qij). Each localized eigenfunction is linearized about a similar spike scaled by the hj0. We thus look for a global eigenfunction composed of similar functions localized in a neighbourhood of each X{. To that end we write, ( c 0 \ 0. (3.63) If we choose the vector c(where c is the vector of scaling coefficients defined in (3.63)) to be any eigenvector of the matrix h 7 - 1 0 h 7 ( m - 1 ) ~ ' s with eigenvalue p then (3.62) becomes, mqp C[(pyy-(f)+ pvPc(j) + Dh c L u^Qdy ) = cA0. (3.64) Thus, we are left with a single scalar equation. Again the form of this equation is similar to (2.19). Using the analysis of this equation from the previous chapter we find that the principal eigenvalue of this system is zero when, 2(s + l)23p 1 Dh > 2' (3.65) We note that each eigenvector of the matrix h 7 1 © h 7 ( m ^ s will give rise to a different global eigenfunction. Thus for each eigenfunction found from (2.24) we may have up to 57 Chapter 3. Finite Inhibitor Diffusion Coefficient n different global eigenfunctions. It will be the smallest eigenvalue denoted by pmin, of h 7 - 10h 7(m - 1)~' s, that will determine the stability of any particular spike configuration. Specifically if r = 2 ( S + l ) 2 ^ t e > l j Dh. 2 then the principal eigenvalue of (3.62) will be zero and the spike configuration will be stable. To illustrate this result we performed some numerical computations with the parameter set (p,q,m,s) = (2,1,2,0). In figure 3.7 and figure 3.8 we plot T versus Dh for a 2 spike and a 3 spike configuration in which the spikes are positioned at (—0.5,0.5) and (-2/3,0,2/3), respectively. The point at which V crosses 1/2 determines the critical value of Dh at which the stability of the spike configuration changes. Figure 3.7 predicts r Figure 3.7: 2 spikes. that a configuration of 2 spikes centered at —0.5 and 0.5 will be stable for Dh < .57 and 58 Chapter 3. Finite Inhibitor Diffusion Coefficient o i I I I I I I I I I I l 0 0 . 1 0 . 2 0 . 3 0 . 4 _ 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1 Figure 3.8: 3 spikes. unstable for larger values of Dh- Figure 3.8 predicts that a configuration of 3 spikes centered —2/3, 0 and 2/3 will be stable for values of Dh < .18 and unstable for larger values. To check these results, we performed numerical simulations on the full system (1.7) using P D E C O L . These computations showed that the two and three spike system are stable for values of Dh < 0.33 and DH < 0.13 respectively. This discrepancy may be due to difficulties in running a simulation near a bifurcation point. The results found for the case of Dh —> oo may be recovered by examining l imrj^oo 0 . We begin by examining limrjJl->0o A. We write r)oti as a power series in D^1, Vo,i = V{o}+D^+0(Df). (3.67) Using the Taylor series expansion of the coefficients of the matrix A, we get the following 59 Chapter 3. Finite Inhibitor Diffusion Coefficient first order system, / i i KAXI 1 KAX2 KAX2 0 Thus f$> K A I I 1 + 1 KAXS K,AX3 1 + KAxn-2 K,Axn-2 K A X „ _ 1 ^ K.Axn-l K A X „ _ I 1 fcAx„ ( do) ^ Vo,0 \ 7 7 ( 0 ) / = f)n-io = ^o°^ To determine fj^ we need to look at the next order system, I KAXI 1 V 0 where, b0 = K2 h = K2 6 n _ x = K2 KAXI '/o,o l K,AX2 KAX2 0 + KAXS KAXS K A I „ - 2 « A l „ - 2 K A X „ _ I + 1 0 1 nAxn-KAX„-KAX„-I / \ bn-l J (3.68) (Axi + Axi+1) - ^(AXi + Axi+1) - sh]^s-l2f3yj ^ + m / ^ ™ " 1 ^ / «rVi,ody I i \ r<x> r J —oo /oo ur'kody -oo K We note that the coefficient matrix of (3.68) has a determinant of zero and a null space spanned by the vector (1 , . . . , 1) T . Thus for a solution to this problem to exist, we require that the left hand side vector is orthogonal to the null space. Taking the dot product of 60 Chapter 3. Finite Inhibitor Diffusion Coefficient (1 , . . . , 1) T with (b0, • • • , 6 n - i ) T and setting the result to zero gives (using the fact that £r=~i l A ^ = 2), "~1 poo J —oo ~ ( 0 ) Vo = 1 2^ + 2snPhlmEs- i = Q (3.69) Therefore © must be tending to the matrix, / l 2K2 + 2sfi \ (3.70) \ 1 ••• 1 / as 7J/i —>• oo. The matrix (3.70) has one eigenvalue of n, with multiplicity 1 and corre-sponding eigenvector (1, • • • , 1) T . The other eigenvalue is 0 with a multiplicity n — 1. The corresponding eigenvectors are of the form e^-, where has a -1 in the ith posi-tion, a 1 in the jth position and zeros elsewhere. The zero eigenvalue implies that the related eigenfunctions will have the same eigenvalues as the local operator. These forms of eigenfunctions correspond to those discussed in section (2.6). 61 Chapter 4 A Spike in a Multi-Dimensional Domain The preceding analysis was carried out in one dimension only. As a first attempt to treat the multi-dimensional case, we will analyze the slow motion of a one-spike solution to the Gierer Meinhardt system in the weak coupling limit Dh —> oo in a multi-dimensional setting. In this limit, we first construct a quasi-equilibrium solution. We then analyze the stability properties of this solution by examining the spectrum of the eigenvalue problem associated with the linearization about this solution. An exponentially small eigenvalue will be shown to be the principal eigenvalue for this linearization. Finally, we use the projection method to derive an equation of motion for the center of the spike. We remark that since some of the calculations below will parallel those in Chapters 1 and 2 rather closely, some of the analysis below will be covered briefly. The non-dimensionalized Gierer Meinhardt system in a domain Q G M.N is Ap At = e2AA- A + —, in Q (4.1a) Am ThHt = DhAH-{iH+—, in Q, (4.1b) An = 0, Hn = 0 on dQ. (4.1c) Here differentiation with respect to n represents the normal derivative on the boundary. As in Chapter 1, we re-scale this system to ensure that the amplitude of the spike is 0(1) as e —> 0. To derive the correct scaling we follow the analysis leading to (1.12), with 62 Chapter 4. A Spike in a Multi-Dimensional Domain the exception that the integration in (1.12) will be replaced by an integration over an N dimensional domain. Thus, (1.13) is replaced by -vh = -vam + vhs + N. (4.2) Solving (1.11) and (4.2), we get = N q = N { P ~ l ) (A *\ Ua {l-p){l + s) + mq, V h (l-p)(l + s ) + m q { -6) Therefore, upon introducing the new variables a and h by A = e~Vaa and H — e~"hh, we obtain the following scaled Gierer-Meinhardt system, which is analogous to (1.15a): at = e 2 Aa — a + —, in Q, (4.4a) p rht = DhAh-u.h + e _ y v — , in O (4.4b) an = 0, hn = 0 on <9Q. (4.4c) We now examine this system in the weak coupling limit Dh —> oo. Following the analysis of Chapter 1, we expand h as a power series in D^1, h^hv + D^hx + OiDf). (4.5) Substituting (4.5) into (4.4) and collecting powers of D^1 we get the following problems for h0 and hi. Ah0 = 0 in fl, (4.6a) am A/i i = Thhot + fiho - £~NTT I N ^> (4-6B) h0„ = 0 on dQ, (4.6c) hln = 0 on dfl. (4.6d) ^From (4.6a) and (4.6c), we conclude that h0 = h0(t), and so h0 is spatially homogeneous. To determine h0 we apply a solvability condition to (4.6b) and (4.6d) to get the following 63 Chapter 4. A Spike in a Multi-Dimensional Domain ordinary differential equation for h0(t): rhh0 + pho r am / = 0- (4-7) Jn no \Cl\ Here h0 = dh0/dt and |f2| is the volume of f2. As was the case for the one-dimensional system, we expect that the dynamics of h are much faster than those of a. Hence, we set h0 = 0 in (4.7) and solve for ho to get i This gives us the Shadow System for the Gierer Meinhardt equations in the weak coupling limit, ap at = e2Aa - a + in Q, (4.9a) ho ho = [ amdx) , (4.9b) an = 0 on an. (4.9c) We now construct a quasi-equilibrium solution aE for (4.9). This is done in a similar manner as in the one-dimensional case, except that the quasi-equilibrium solution will be radially symmetric about the center of the spike. Thus, we look for a solution to (4.9) in all of in the form a = aE(x) = K0yuc(p), p = e _ 1 | x - x 0 | , 7 = q/(p - 1). (4.10) The function u c(p), called the canonical spike solution, is radially symmetric about the origin and it decays exponentially as p —> 00. It satisfies ul + ^^u'c + uc-upc = Q, (4.11a) <(0) = 0 and uc(oo) = 0, (4.11b) uc(p) ~ ap{1-N)/2e-p, as p -»• 00 . (4.11c) 64 Chapter 4. A Spike in a Multi-Dimensional Domain In terms of this solution, the quasi-equilibrium solution is given by aE(x) = hluc (e x|x - x0|) , (4.12a) 0 " U N Jo h ° = 777% / < V V - i d p • (4.12b) Here QN is the surface area of the unit iV dimensional sphere. Recall that in the one-dimensional case and with p = 2 we have the exact solution of uc(p) = | sech 2 ( | ) . To find numerical solutions for uc(p) in other dimensions, we will treat N as a real parameter, and use iV (and p for p ^ 2) as continuation parameters. We can use the far field asymptotic behavior (4.11c) to obtain the boundary condition u'c = uc, which we impose at some large value p = pL in our numerical computations of (4.11). The computations are done using C O L N E W . In Fig. 4.1 we plot the numerically computed solutions uc(p) for N = 1,2,3 when p = 2. Chapter 4. A Spike in a Multi-Dimensional Domain Again, we note that aE will satisfy the steady-state problem for (4.9a), but will fail to satisfy the no flux boundary condition (4.9c) by only exponentially small terms for any value of x 0 G f2. Thus, we expect that the spectrum of the eigenvalue problem associated with the linearization about aE contains exponentially small eigenvalues. We now obtain the eigenvalue problem for this linearization by introducing <j> and 77 defined by o(x, t) = a £ (x) + e A V ( x ) , (4.13a) /i(x, t) = h0 + eA*77(x). (4.13b) Here (j) <C aE and 77 <C h0. We substitute (4.13) into (4.4) to arrive, after a lengthy calculation, at the following implicit eigenvalue problem, pNUN(p(s + 1)4- AT) JN (4.14a) (j)n = 0 on dfl. (4.14b) Here, j3^ is defined by POO 0N= / u™-lpN~ldp. (4.15) Jo In (4.14), uc = u c [e _ 1 |x — x 0 | ] . Thus, we will only seek eigenfunctions that are localized near x = x0. These eigenfunctions are of the form 0(y) = 0 ( x o + ey), y = e- 1(x-x 0). (4.16) Therefore, we can replace fl by W1 in (4.14a) and impose a decay condition for 4> as |y| —¥ 00. This gives us the eigenvalue problem for the infinite domain L£4> = Ay4> + (-1 + pup-l)i> - ~—^L_— / unc-14>dy = \4>, in RN, PNilN(p(s + 1) + Ar) JRN (4.17a) 4>-¥0 as | y | - ) - co . (4.17b) 66 Chapter 4. A Spike in a Multi-Dimensional Domain In this problem uc = i t c ( | y | ) . If, in addition, we consider an eigenfunction that is radially symmetric (i. e. 0 = fap), where p = |y|), then (4.17) reduces to (4.18a) </>->• 0 as p - » o o , (4.18b) where A p 0 = fa + (N — l)p~lfa Since r is typically very small, we will use the simplifying approximation that r = 0 for the remainder of the analysis in this section. We first note that the function fa = dyiuc(\y\) for i = 1,..,N satisfies (4.17). Here Hi is the i ^ n coordinate of y. This follows from the combined effects of translation invariance and the vanishing of the integral in (4.17) by symmetry considerations. Thus, this problem has a zero eigenvalue of multiplicity N with corresponding eigenfunctions <fi — dyiuc(\y\) for i = 1,.., N. These eigenpairs will be perturbed by only exponentially small terms as a result of the finite domain. Hence, there are N eigenvalues of (4.14) that are exponentially small. The goal is to determine whether these are the principal eigenvalues of (4.14). If the" non-local term is absent in (4.14) then, since each of these eigenfunctions has one nodal line, we know that these eigenvalues are not the principal eigenvalues. When the non-local term is absent, there is exactly one principal eigenvalue for (4.14) and this eigenvalue is positive and bounded away from zero. Thus, no meta-stable behavior can occur when the non-local term in (4.14) is absent. By introducing the effect of the non-local term in a gradual way, we will show that this positive eigenvalue will cross through zero and become negative. Hence, the effect of the non-local term will be to ensure that the exponentially small eigenvalues are the principal eigenvalues for the non-local eigenvalue problem (4.14). We also note that When the non-local term is absent in (4.14), the principal eigenfunction is radially symmetric. Thus, we will track this eigenfunction as a function of a continuation parameter that gradually introduces the effect of the non-local term in (4.14). 67 Chapter 4. A Spike in a Multi-Dimensional Domain Therefore, we must first compute the eigenvalues and eigenfunctions of the radially sym-metric problem (4.18) (with r = 0). This calculation is done by repeating the procedure used in one dimension, in which the non-local behavior is introduced gradually through a continuation parameter 5 Ls4> = AP4> + (-l+ pvP-l)j> - ^ 7 M ^ <~XhN-X dp = Xfa p > 0, (4.19a) 0-» 0 as p-+oo. (4.19b) We need to determine the eigenvalues of this problem as a function of 5 and to confirm that its first eigenvalue X0(5) has a negative real part when 5 = 1. To solve this eigenvalue problem numerically we take N as a continuation parameter as well as 5 (we could use pas a continuation parameter for p ^ 2). In Fig. 4.2 and Fig. 4.3 we plot the first two eigenvalues Xo(5) and A2(5) of (4.19) as a function of 5 for iV = 2 and N = 3, respectively. These computations were done using C O L N E W . These plots clearly indicate that A0(<5) crosses through 0 before 5 = 1. At some value of 5, A 0 and A 2 collide and become complex. To track the eigenvalues past the point where they become complex, we use the same technique as in the one-dimensional case. The differential operator is approximated by a matrix and the eigenvalues of the matrix are then approximations of the eigenvalues of the differential operator. Using this numerical procedure, we give numerical values for the real and imaginary part of XQ(5) in Table 4.1. This table shows that the real part of A 0 is negative when 5 = 1. 4.1 An Exponentially Small Eigenvalue We will now use a boundary layer analysis to construct a composite approximation to the eigenfunctions corresponding to the exponentially small eigenvalues of (4.14). The corresponding eigenfunction is well approximated by dXiuc in the interior of the domain and has a boundary layer correction term near dfl in order to satisfy the no-flux boundary condition on dfl. In order to resolve the boundary layer we must define a local coordinate 68 Chapter 4. A Spike in a Multi-Dimensional Domain Figure 4.2: \o(S) and A2(5) versus 5 in E 2 for the parameter set(p = 2, q = 1, m = 2, s = 0). system. Let 77 represent the distance from a point in Q, to <9Q, where 77 < 0 corresponds to the interior of Q,. Let £ correspond to the other N — 1 orthogonal coordinates. To localize the region near dQ, we let 77 = e~lf). The eigenfunction on the finite domain can then be approximated by, <j)i = Ci (dXiuc (e x|x - x0|) + fc) , (4.20) where Ci is a normalization constant and fc is a boundary layer correction term. Using the fact that uc is exponentially small near d£l we get the following boundary layer 69 Chapter 4. A Spike in a Multi-Dimensional Domain Figure 4.3: X0(8) and X2(8) versus S in R3 for the parameter set (p = 2, q = 1, m = 2, s 0). problem dmfa - fa = 0, n < 0, df) dvfa = - ^ ( 5 S i M c ) | ^ = 0 — , on ?7 = 0, a function of C 0i —> 0 as 77 —> — 0 0 . (4.21) (4.22) (4.23) We require that fa —> 0 as r\ —> — 0 0 to match to the outer solution. Define </i(C) to be the right side of (4.22), ft(C) = -^ (a s i u c ) | , = o . (4.24) Then, the solution for fa is fa = 9i(C)ev. (4.25) 70 Chapter 4. A Spike in a Multi-Dimensional Domain Thus, the composite asymptotic solution for the eigenfunction is <t>i = Ci[dXiuc + gi{Qe1f\, i = l...,N. (4.26) In order to complete our asymptotic estimate of the exponentially small eigenvalues, we apply Green's identity to fc and dXiuc to get the following relationship: Xi(dXiuc, fc) = - e 2 / (t>idn{dXiuc)dS + (L*dXiuc, fc). (4.27) Jan Here L*£v EE e2Av - v + ul~lv - J <v dx . (4.28) We will now estimate each term in (4.27). Since dXiuc is an exact solution to the local problem, we have that, K{dXluc) = - ^^1 j f ^ t i c d x . (4.29) Next, since uc is radially symmetric and localized to a small region in the interior of Q,, it is clear that fnupdXiucdx. = §uupdx.ucdx., = 1... N. Thus, we may write the expression above as, ,m—1 " N l_c 2N3(s +1) A n application of the Divergence Theorem results in, .m— 1 L_c 2N{3(s + On the boundary of £2, uc [e~x|x - x 0 | ] ~ a 6(JV-i)/2 r(i-JV)/2 e- £-i|x-x 0|_ Therefore, the integral in (4.31) will be exponentially small. If the boundary of Q, is smooth, we can estimate the integral to get the following bound: \r*(8 7/ )\ < \ m \ p+i (JV-IMP+D/2 (I-AQ(P+I)/2 6 - i ( p + i ) P n „ 2 ) I M < W I < 2N(3{s + 1){p+1f Po (4.3^ j 71 K (9Xiuc) = - . ^ r ^ / E«Sa X iu cdx. (4.30) J n i=i L : ( 0 ^ ) = - ^ _ J ^ ) ^ (4.31) Chapter 4. A Spike in a Multi-Dimensional Domain Here p 0 = dist(x 0, dtt). Therefore, we have Thus, X £ ( » - i ) ( p + i ) / 2 ^ - « ) ( P + i ) / 2 ( , - ( P + i ) ( w - ' / 2 [ u?-iaXiucdx. (4.33) Jn A similar procedure shows that mq \dQ,\2 \(L:dXi,fa)\< 2N2f3(s + l) (p+l)m x a ( p + l ) m £ ( l - i V ) ( p + l ) m / 2 ^ ( l - J V ) ( p + l ) m / 2 e - ( p + i ) ( m - i ) p 0 / ^ '±34) Thus, this quantity is exponentially small. We will show that it is exponentially smaller than the other terms of (4.27). Therefore, we can ignore it. To proceed with the analysis, we need estimates for the eigenfunctions on the boundary. The first step is to find gj(() in (4.24). Let x0j represent the jth component of x0 and let Xj = Xj(Q be a parameterization of the boundary. So, setting r = \x — x0\, we apply the chain rule, which gives 9 i ~ - { X i ^ K ( r A ) r - n ) , (4.35) where n is the outward unit normal to £1. Since uc(p) ~ ap^~N^2e~p as p —> oo we get that, Xj XQJ ) r - ( 1 + J V ) / 2 e - r / e r • n, on 5 0 . (4.36) We can now combine (4.20) with (4.36) to get an asymptotic approximation of fa on the boundary. In this way, we find fa~ - C i e ( J V - 3 ^ 2 a ( o ; i - r r i o ) r - ( 1 + i V ) / 2 e - r / e ( l + r-ri), on OQ. (4.37) This expression will be used in the integral term in (4.27). Now we estimate the left hand side of (4.27). Since fa and dXiuc are exponentially small outside of a neighbourhood of x = xo, this inner product will be dominated by these 72 Chapter 4. A Spike in a Multi-Dimensional Domain localized functions. Using a Laplace-type approximation, we can approximate the inner product to get (0 X i « c > fa) ~ % j(Vc(r/e))2 ( x - i ^ i \ d x = 9^1 j u'c{p)2pN-1 dpd9, (4.38) 6 Jn \ r J i V Jw where 6 represents the N — 1 angular co-ordinates. Since the integrand is independent of 6, we can define /3N = f0°° u'c(p)2pN~1dp, to simplify (4.38) to (dXiuc, fa) ~ CjeN-2£lNpN/N. (4.39) Here fl]^ is the surface area of the n-dimensional unit sphere. We may now find Cj by using the normalization relation fn fa dx = 1 to obtain, N 1/2 d = ( ^ — ) e^'2. (4.40) \PNSIN) Finally, we get our asymptotic estimate of Ai by substituting (4.39), (4.37) into (4.27) and using the estimate dn(dXiuc) ~ ae^ J V _ 5 ^ 2 r~^ i V + 1 ^ 2 e _ r / e , on dQ. In this way, we get ^ ~ S^- [ fa - x 0 i ) 2 r - ( 1 + A V 2 r / e ( r • n)(l + r • n) dS. (4.41) PN^U Jon As a consistency check we observe by comparing the asymptotic orders of the two terms on the right side of (4.27) that the second term is asymptotically negligible compared to the first term. This surface integral may be further simplified by using a multi-dimensional Laplace technique. We let x m be the point on dQ, where rm = dist(xo, dfl) is minimized. Assume that x m is unique. If we parameterize the boundary near Cm (where x(Cm) = x m ) such that each Q corresponds to arclength along one of the principal directions through £ m , then for any smooth F(r), we have(see [15]), „ / v (N-l)/2 / r1-NF(r)e-2r^dS= — F(rm)H(rm)e-2r^e, (4.42) Jan \rmJ 73 Chapter 4. A Spike in a Multi-Dimensional Domain where • H(rm) = (1 - rm/Ri)~1/2(l ~ rm/R2)-1/2 • • • (1 - r j R ^ y 1 ! 2 . (4.43) Here Rj are the principal radii of curvature of dfl at x m . This result assumes the non-degeneracy condition Rj > rm, j = 1,... , N — 1 holds. Using equation (4.42) we have the following asymptotic estimate for the exponentially small eigenvalue, A " - f ^ ( r ) M / 2 ( ^ ) 2 f f ( - ) e " 2 ' " / 2 - < 4 4 4 > PjVliJV \'m/ \ 'm / We have used the fact that at xm, r will be in the same direction as the normal vector and thus r • n = 1 at this point. Now we may examine the dynamics of the iV-dimensional spike system. We linearize about a moving spike by writing, a(x, t) = KQUC ( e - 1 | x - x 0(£)|) + w(x,t), (4.45) where w < h^uc (e _ 1 | x - x 0(t)|) and wt < e~lr~lhlu'c (e _ 1 | x - x 0 ( i) | ) 2~2?=i x0i(xi-x0i). Substituting (4.45) into (4.9), we obtain, N Lew = e~1r~1h'^u'c(e~1\x-ico{t)\)'^2xoi(xi - x0i), in fi, (4.46a) i=l 9 n w = - £ - 1 / i X ( € ~ 1 | x - x o ( t)|)nT, on dQ. (4.46b) Now we wish to isolate the behaviour caused by the exponentially small eigenvalues. We will expand w as an eigenfunction expansion. As in chapter 1, since our operator is not self-adjoint, we will use the eigenfunctions of the local operator. Again we will utilize the fact that the eigenpairs of interest are common to the local and adjoint operators, as well as L£. We will refer to the eigenpairs of the local operator and the adjoint operator by (Aj, fc) and (A^, </>*), respectively. We now expand w in terms of the eigenfunctions of the local operator as they form a complete orthonormal set, ^ E ^ " , (4-47) »=i ' 74 Chapter 4. A Spike in a Multi-Dimensional Domain The coefficient D ; in (4.47) is given by Di(t) = Xi(w, fa) e2 / dnwfadS + e lr lhl\u'c(e *|x - x0(<)|) y^jxoj(xj - x0j),fa } Jan V j=i J = - e 2 / dnw fa dS + hi £ x0jCi / dXjuc (e _ 1 | x - x 0 ( i) | ) <9Xiuc ( e - 1 | x - x0(<)|) dx. 7an - = 1 (4.48) Clearly (<&, </>j) = 0, for i ^ since uc is a radially symmetric function, and fa ~ 5 Z i « c -Then (4.48) simplifies to Di(t) = - e 2 / dnwfadS + h1x0iCi / [d X i « c (e _ 1 | x - x 0 ( i ) | ) ] 2 dx, Van i n 7 . (4.49) = - 6 2 / dnwfadS + ^ Jan Since A; is exponentially small, we may apply the solvability condition that Di(t) = 0 , for i = 1,... N. This will result in the following N coupled ordinary differential equations for the coordinates of the spike location x 0 : i <M = e2S f dnwfadS, i = l,...,N. (4.50) "0 Jan We can evaluate the right hand side in the expression using (4.46b) and our estimates for fa in (4.37) and for uc in the far field. This gives, xoi ~ eN-lC2a2 [ (xi - x0i)r-Ne-2r/e(l + r • n)r • n dS Jan ^feNo_\ f (x._Xo.y-Ne-2r/e(1 + T . n ^ . n d S V PN^N ' Jan (4.51) The equation above tells us that the spike will be repelled by the boundary and thus will tend towards the center of Q, (see [15]). This is the main result of this section. We can obtain the equilibrium position of the spike by setting xo = 0 and solving the resulting transcendental equations. A geometrical interpretation of this result shows that 75 Chapter 4. A Spike in a Multi-Dimensional Domain the center of the spike will tend towards an 0(e) distance of the center of the largest inscribed sphere in Cl. This result requires Cl to be convex (see [15]). We can simplify the above equations by applying equation (4.42), which results in the following equations of motion, V P / v r i i j V / \ ' m / I m One must be careful when applying the above equation, for as the spike moves the value of x m may change. 76 Chapter 4. A Spike in a Multi-Dimensional Domain s Ao Ao in E 2 in E 3 0.00000 1.6388 2.3703 0.05000 1.4814 2.1588 0.10000 1.3231 1.9456 0.15000 1.1638 1.7304 0.20000 1.0030 1.5125 0.25000 0.84032 1.2910 0.30000 0.67516 1.0646 0.35000 0.50641 0.83098 0.40000 0.33218 0.58554 0.45000 0.14857 0.31741 0.50000 -.055026 -.019898 0.55000 -.37526 -.33843 + 0.29744? 0.60000 -.48239 + 0.24569? -.44368 + 0.45028? 0.65000 -.56115 + 0.33165?; -.54978 + 0.54508? 0.70000 -.64059 + 0.38475? -.65696 + 0.60964? 0.75000 -.72097 + 0.41770? -.76550 + 0.65310? 0.80000 -.80268+ 0.43510i -.87584 + 0.67970? 0.85000 -.88640 + 0.43886? -.98857 + 0.69170? 0.90000 -.97333 + 0.42959? -1.1045 + 0.69037? 0.95000 -.10657 + 0.40726?' -1.2249 + 0.67652? 0.10000 -.11678 + 0.37248? -1.3513 + 0.65089? Table 4.1: S and A 0 in E 2 and E 3 for the case of (p = 2, q = 1, m = 2, s = 0). 77 Chapter 5 Conclusions Previous studies have shown that scalar reaction diffusion systems can not support stable spike-like solutions. However, two species reaction diffusion systems have been shown to exhibit stable spike like solutions. To examine the cause of this stability change, we linearized about a one-spike solution and studied the resulting linear operator. The difference between the linear operator resulting from a scalar reaction diffusion equation and a two species system is the addition of a non-local term. It is this integral term that allows the system to be stabilized. When we consider n-spike solutions, the situation is more complex. One might assume that the stability of an isolated spike of an n spike solution will only depend on the local levels of inhibitor. This has proven not to be the case. For if it were, in the limit of Dh —>• oo, an n-spike solution with too many spikes would tend to zero globally in space as each spike would feel the same level of inhibition. However, it is observed numerically that some spikes persist while other vanish. Thus, the determination of the stability of n-spike solutions is very different for the two cases studied. The stability and dynamics of the Gierer Meinhardt equations were studied first for the limiting case of Dh oo and e —> 0. For this case it was found that one-spike solutions are stable, as the equilibrium level of inhibitor for one spike is sufficiently small to allow the spike to persist. This fact is reflected in the spectrum of the linearized operator resulting from a linearization about a one-spike solution. However, in this limit it is 78 Chapter 5. Conclusions impossible to find n-spike solutions, with n > 1, which are stable. This is reflected by the nature of the eigenfunctions resulting from a linearization about an n-spike solution. For this limit one can always find an eigenfunction with two localized extrema that has an eigenvalue bounded above from zero. As the modes corresponding to these positive eigenvalues grow, some spikes will vanish, while others will persist. It is conjectured that the winner of this spike competition will be determined by the positioning of the spikes. I suspect that the determination of this winner will be exponentially sensitive. The projection method was used to quantify the motion of a one-spike solution. The presence of an exponentially small principal eigenvalue is used to imply a limiting solv-ability condition on the system, which in turn provides an equation of motion for the center of the spike. However, the full numerical results appear to move on a much faster time scale than predicted by the asymptotics. Difficulties in simulating an exponentially sensitive system may be responsible for this discrepancy. In the case Dh = 0(1) as e —>• 0 it is possible to have stable n-spike solutions. In this case the heights of the spikes will be dependent on their locations. A quasi-equilibrium solution may be found by considering each spike as a point source of inhibition. The heights of the spikes and the solution to the inhibition profile may then be computed simultaneously. The stability is dependent on both the value of Dh and the positions of the spikes. In this limit, the stability of an n-spike solution will depend on the eigenvalues of a matrix that couples the localized eigenvalue equations. The eigenvectors of this matrix form scalings of the local eigenfunctions, which when added together form a global eigenfunction. The results from numerical simulations of the full system agree well with the asymptotic stability predictions. There is a slight discrepancy, but this is likely due to problems in quantifying the stability of the system using full numerics. Qualitatively, the numerical and the asymptotic solutions agree both in motion and sta-79 Chapter 5. Conclusions bility. Since both of these phenomena are controlled by small eigenvalues, the numerical problem will be ill-conditioned and accurate results are difficult to obtain. Numerically simulated unstable solutions, which are the result of an eigenvalue of order one, agree very well with the predictions (see Table 2.3). There are many questions that remain to be answered about the Gierer Meinhardt system. For the case of large Dh, the question of determining which spike of an n-spike solution will persist is still to be answered. I believe that the non-local effects of the operator will have to be considered to determine which eigenfunction will have the largest eigenvalue. For finite values of Dh, n-spike solutions can be stable. It should then be possible to find the equations of motion governing an n-spike system. There are few results for the Gierer Meinhardt system in a multi-dimensional setting. Only the case of one spike with large Dh has been analyzed so far. The extension of the case to finite Dh and to multiple spatial dimensions should, conceptually, present no problems. However, the technical details of incorporating an arbitrary boundary will require more careful consideration. 80 Bibliography Bruce Alberts, Dennis Bray, Julian Lewis, Martin Raff, Keith Roberts, and James D. Watson. Molecular Biology of the Cell. Garland Publishing, Inc, third edition, 1994. Uri M . Ascher, Steven J. Ruuth, and Brian T.R. Wetton. Implicit-explicit methods for time-dependent partial differential equations. SIAM J. Numer. Anal., 32(3):797-823, June 1995. Pedro Freitas. A nonlocal Sturm-Liouville eigenvalue problem. Proc. Roy. Soc. Edinburgh, 12 A: 169-188, 1994. Jr. G. L. Lamb. Elements of Soliton Theory. Wiley Interscience, 1980. A. Gierer and H. Meinhardt. A theory of biological pattern formation. Kybernetik, 12:30-39, 1972. David M . Holloway. Reaction-diffusion theory of localized structures with application to vertebrate organogenesis. PhD thesis, University of British Columbia, March 1995. James P. Keener. Activators and inhibitors in patern formation. Stud. Appl. Math, 59:1-23, 1978. Wei-Ming Ni and Izumi Takagi. On the shape of least-energy solutions to a semilinear Neumann problem. Comm. Pure Appl. Math, XLIV:819-851, 1991. Wei-Ming Ni and Izumi Takagi. Point-condensation generated by a reaction-diffusion system in axially symmetric domains. Japan J. Indust. Appl. Math, 12(2):327-365, 1995. Y . Nishiura. Global structure of bifurcating solutions of some reaction-diffustion systems. SIAM J. Math. Anal., 13(4):555-593, 1982. Steven J. Ruuth. Implicit-explicit methods for reaction-diffusion problems in pattern formation. J. Math. Biol, 34:148-176, 1995. I. Takagi. Point-condensation for a reaction-diffusion system. J. Differential Equa-tions, 61:208-249, 1986. A . Turing. The chemical basis of morphogenesis. Phil. Trans. Roy. Soc. B, 32, 1952. Michael J. Ward. Eliminating indeterminacy in singularly perturbed boundary value problems with translation invariant potentials. Stud. Appl. Math., 87(2):95-134, 1992. 81 Bibliography [15] Michael J . Ward. An asymptotic analysis if localized solutions for some reaction-diffusion models in multi-dimensional domains. Stud. Appl. Math, 97(2): 103-126, 1996. [16] Michael J. Ward. Metastable bubble solutions for the Allen - Cahn equation with mass conservation. SIAM J. Appl. Math., 5(1247-1279), 1996. [17] Michael J. Ward. Boundaries, Interfaces and Transitions, chapter Dynamical metastability and singular perturbations. American Math Society, 1997,to appear. 82
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Meta-stability of the Gierer Meinhardt equations
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Meta-stability of the Gierer Meinhardt equations Iron, David 1997
pdf
Page Metadata
Item Metadata
Title | Meta-stability of the Gierer Meinhardt equations |
Creator |
Iron, David |
Date Issued | 1997 |
Description | A well-known system of partial differential equations, known as the Gierer Meinhardt system, has been used to model cellular differentiation and morphogenesis. The system is of reaction-diffusion type and involves the determination of an activator and an inhibitor concentration field. Long-lived isolated spike solutions for the activator model the localized concentration profile that is responsible for cellular differentiation. In a biological context, the Gierer Meinhardt system has been used to model such events as head determination in the hydra and heart formation in axolotl. This thesis involves a careful numerical and asymptotic analysis of this system in one dimension for a specific parameter set and a limited analysis of this system in a multidimensional setting. Numerical analysis has revealed that once the spikes form they continue to move on an extremely slow time scale. This type of phenomenon is a general indicator of meta-stable behaviour. By perturbing off of an isolated spike solution an exponentially small eigenvalue of the linearized operator was found. This small eigenvalue accounted for the extremely slow motion found numerically and thus was used to obtain an equation of motion for the location of the spike. The Gierer Meinhardt system is analyzed in the limit of small activator diffusivity for both a finite inhibitor diffusivity and for an asymptotically large inhibitor diffusivity. In this thesis, the mathematical techniques used include the method of matched asymptotic expansions, spectral theory and numerical computations. |
Extent | 3344968 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-24 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080011 |
URI | http://hdl.handle.net/2429/6402 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1997-0419.pdf [ 3.19MB ]
- Metadata
- JSON: 831-1.0080011.json
- JSON-LD: 831-1.0080011-ld.json
- RDF/XML (Pretty): 831-1.0080011-rdf.xml
- RDF/JSON: 831-1.0080011-rdf.json
- Turtle: 831-1.0080011-turtle.txt
- N-Triples: 831-1.0080011-rdf-ntriples.txt
- Original Record: 831-1.0080011-source.json
- Full Text
- 831-1.0080011-fulltext.txt
- Citation
- 831-1.0080011.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080011/manifest