A Computer-Assisted Study of Red Coral Population Dynamics
Sayomi Kamimoto, Hye Kyung Kim, Evelyn Sander, Thomas Wanner
AA Computer-Assisted Study of Red Coral Population Dynamics
Sayomi Kamimoto , Hye Kyung Kim , Evelyn Sander , Thomas Wanner
1. Department of Mathematical Sciences, George Mason University,Fairfax, VA 22030, USA2. School of Mathematics, University of Minnesota,Minneapolis, MN 55455, USASeptember 30, 2020
Abstract
We consider a 13-dimensional age-structured discrete red coral population model vary-ing with respect to a fitness parameter. Our numerical results give a bifurcation diagramof both equilibria and stable invariant curves of orbits. We observe that not only for lowlevels of fitness, but also for high levels of fitness, populations are extremely vulnerable, inthat they spend long time periods near extinction. We then use computer-assisted proofstechniques to rigorously validate the set of regular and bifurcation fixed points that havebeen found numerically.
AMS subject classifications:
Primary: 37G15, 37M20, 65G20, 65P30; Secondary:37B35, 37C70, 65G30, 92D25, 92D40.
Keywords:
Bifurcations, Computer-Assisted Proofs, Red Coral, Age-Structured Popu-lation Models, Interval Arithmetic, Rigorous Validation
Contents a r X i v : . [ m a t h . D S ] S e p Introduction
Coral plays an important role in the marine ecosystem, and coral reefs provide habitatsto many sea animals and protect coastlines from breaking waves and storms. Red coralis a long-lived, slow-growing species, dwelling on Mediterranean rocky bottoms. Red coralpopulations are at risk due to both global climate change and overharvesting [4]. Bramanti,Iannelli, and Santangelo [3, 17] investigated red coral populations by scraping samples offthe coast of Italy in Calafuria in the Western Ligurian Sea (43 ◦ (cid:48) N, 10 ◦ (cid:48) E, Italy, at adepth between 20 and 45m depth) and observing their growth rate over a four-year period.They used this data to construct a Leslie-Lewis transition matrix, a static life table, anda 13-dimensional dynamical population model. Using this model, they studied populationtrends by comparing small young colonies and bigger older colonies. However, they onlyconsidered a small range of population trends. In the current paper, we present a systematicstudy of this coral population model, shedding light on the long-term dynamics of the redcoral populations. We can see the long-term effect of change in reproduction fitness. Weestablish the equilibrium structure and bifurcation points for the model, find a set of stableperiodic invariant cycles, and show that for a large range of reproduction fitness these cyclesget close to population extinction.In addition to these observations, we present and implement methods which allow usto rigorously validate the model’s equilibrium and bifurcation structure, including both asaddle-node and a Neimark-Sacker bifurcation. These validations use a modification of theNewton-Kantorovitch type method developed in [15, 19, 20]. While the previous versionof this method merely used natural continuation, this paper contains an extension of theseresults in which we consider rigorous validation using pseudo-arclength continuation [9, 10].In addition, we use computer-assisted proof methods to prove the existence of saddle-node andNeimark-Sacker bifurcation points on the equilibrium branch. These methods significantlyextend the range of applications of the constructive implicit function theorem which wasintroduced in [15]. While for the purposes of this paper we restrict ourselves to the case offinite-dimensional Euclidean spaces, the results can easily be adapted to the general Banachspace setting, with little change. Thus, the pseudo-arclength results can be used for examplein the setting of partial differenial equations, such as the setting described in [16]. In otherwords, the present paper presents a functional analytic foundation for using pseudo-arclengthcontinuation in the context of computer-assisted proofs based on the constructive implicitfunction theorem presented in [15].The remainder of this paper is organized as follows. We introduce the age-based red coralmodel in Section 2. In addition, we present a bifurcation diagram of fixed points and stabilityof the model, along with a detailed discussion of oscillations. These results show how even athigh fitness levels, the oscillations lead to extreme vulnerability of the population. Section 3contains a functional-analytic approach to the rigorous validation of the regular branches inthe bifurcation diagram, which is based on a constructive version of the implicit functiontheorem. Subsequently, Section 4 details the validation for the three bifurcation points onthe main fixed point branch; namely, the saddle-node bifurcation in 4.2, the Neimark-Sackerbifurcation in 4.1, and the transcritical bifurction in 4.3. Section 5 contains conclusions andfuture work. 2igure 1: Photographs of red coral colonies. The individual polyps are visible particularly inthe right-hand image. Photos from [1, 8].
In this section we present the red coral population model of Bramanti, Iannelli, and Santan-gelo [3, 17], based on their experimental and field data and a Leslie-Lewis transition matrix.In addition, we describe the dynamics of the model in terms of its bifurcation structure anddiscuss its implications.
A coral population is a self-seeding independent group consisting of polyps , tiny soft-bodiedorganisms related to jellyfish. Polyps form into colonies , which are distinct clusters withpolyps residing on a surface, as shown in Figure 1. A polyp is born to a parent colony in afree-swimming larval stage. At the end of the larval stage, the polyp permanently attachesitself to a colony and cannot move again. The age of a colony has implications in terms ofits size and polyp density. As a result, colony age determines the polyp attachment rate, thelarval birth rate, and the polyp survival rate. Based on these factors, larvae will attach eitherto an existing colony or, especially if there is a high polyp density, recruitment will occur.That is, larvae do not attach to existing colonies, but instead form new colonies. Red coralpolyps can reproduce larvae starting two years after their birth, implying that there is nobirth in a colony less than two years old, since none of the polyps are old enough to reproduce.Reproduction occurs at a discrete time in summer, implying that a discrete population modelis a natural modeling assumption.Based on the setting above, rather than modeling the total large number of polyps in acoral population, the age-based model is a discrete time model for ( x , x , . . . , x d ), where x k is the number of colonies of age group k . The value d is the oldest colony in the population.While in principle this d could be large, in the observations made there was no colony of agegroup greater than 13. The value of x k changes with respect to time (in years), where x nk denotes the number of colonies of age group k at year n . The colony life cycle is displayedin the schematic diagram shown in Figure 2. The downward arrows in Figure 2 indicatethat x k +1 , the number of colonies in age group k + 1, is determined exclusively by the3igure 2: Life cycle of coral population Class k Survival rate S k Fertility F k k in the previous year. This relation is linear with respectto population, with the survival rate constant S k . That is, we have x nk +1 = S k x n − k . Thesurvival rate values are determined by observation, and are given in Table 1, based on [17,Table 2].The upward arrows Figure 2 indicate that recruits may be larvae from any colony of agetwo or greater. Though it is not obvious from the schematic diagram, the recruitment rate isnot linear, and it depends on both the total number of polyps in the colonies, as well as onthe larvae birth rates. Considering that the base variables x k denote the number of coloniesin age group k , the total number of polyps can be deduced from the numbers p k of polypsper colony in a colony of age group k , and the birth rates b k depend on the fertility rates F k given in Table 1. Combined with the observational data in [3], Bramanti et al. have thenderived empirical expressions for the polyp per colony numbers p k and the birth rates b k ,which are given by p k = 1 . k . and b k = F k k . . (1)For our calculations in the present paper, we use these fitting functions rather than theoriginal data, in keeping with the equations in [3]. In addition to the birth rates, the numberof recruits x depends also on a nonlinear function ϕ , which in turn depends on the densityof polyps per unit area. This function ϕ is given by ϕ ( y ) = c e − αy y + c e − βy , with c = 1 . · , c = 1 . · , α = 5 · − , β = 3 . · − , (2)which again is a fit for the observational data in [3]. The shape of this nonlinearity isdepicted in Figure 3. For a small density of polyps, the function ϕ increases with polyp4
500 1000 1500 2000 2500 3000
Polyp density ? Figure 3: The recruits-to-larvae ratio function ϕ plotted with respect to polyp density P .density, whereas too large of a polyp density inhibits the creation of new colonies due tocompetition for resources.We now explain how to compute the polyp population density P . We have already seenthat the numbers p k of polyps per colony in a colony of age group k satisfy the empiricalformulas in (1). Thus, the total number of polyps in age group k is given by p k x k . Nowlet Ω denote the total area of the population site, which was measured to be equal to 36 dm in [3]. Moreover, let x = ( x , x , . . . , x d ) be a column vector giving the number of colonies ofeach age group, and let p = ( p , p , . . . , p d ) denote the vector of polyps per colony in each agegroup. Then the total number of polyps in the (non-recruit!) population Q and the polyppopulation density P satisfy the identities Q = d (cid:88) k =2 p k x k and P = Q Ω . (3)Based on these preliminaries, let x n = ( x n , x n , . . . , x nd ) represent the vector containing thenumber of colonies at year n , and let P be the polyp population density defined in (3). If wenow define L ( λ, x ) = λb ϕ ( P ) λb ϕ ( P ) λb ϕ ( P ) . . . λb d − ϕ ( P ) λb d ϕ ( P ) S . . . S . . . S . . . . . . S d − , (4)5here the bifurcation parameter λ is described below, then our model is given by x n +1 = L ( λ, x n ) x n . (5)The model (4) and (5) is an age-structured, nonlinear, discrete-time dynamical model. Forthe parameter value λ = 1, it is precisely based on the observational data in [3]. Thenonlinearity arises only in the evolution of the variable x , which describes the number ofrecruit colonies. In a slight reformatting of notation, let the function f : R × R d → R d begiven by f ( λ, x ) = L ( λ, x ) x . Then x n +1 = f ( λ, x n ), meaning that the dynamical populationvariation corresponds to the iteration of the parameter-dependent nonlinear map f .We still have to justify the introduction of the bifurcation parameter λ in the above formu-las. Previous work concentrated on the effect of varying the biologically relevant reproductivenumber R , the total number of larvae produced by a single colony during its life span. Thisparameter is directly proportional to λ , as we will show in Section 2.3. The birth rate param-eters b k in the above equation are determined by observation of a specific coral populationover a small time period. In order to consider a population model in which the populationis placed under stress, such as in the case of climate change, it is necessary to change theparameters beyond what has been observed. While we could also consider modification ofother parameters, we choose to follow along the lines of [3] and vary the birth rates, makingthe assumption that every birth rate parameter will be equally affected. Therefore, in oursubsequent analysis, for every k we let the birth rate be given by λb k , a fixed scaling factorcompared to the originally observed birth rate. We now consider the set of fixed points for the coral population model, given by the nonlinearfunction f defined above, and how this set changes as a function of the parameter λ . Thatis, we wish to determine the set of all pairs ( λ, x ) ∈ R × R d such that f ( λ, x ) = x . As it turnsout, this can be reformulated equivalently as a one-dimensional problem. To see this, assumethat we have x = f ( λ, x ). Then for all indices k = 1 , . . . , d − x k +1 = S k x k . Usingthese statements iteratively, one readily obtains x = S x , x = S S x , x = S S S x , . . . x d = S d − · · · S S x . Thus, for all k = 2 , . . . , d we have x k = a k x , where one uses the abbreviation a k = k − (cid:89) i =1 S i , (6)and we further define a = 1 then one also has x = a x . Since we can write each com-ponent x k for k ≥ x alone, the fixed point problem is a one-dimensionalproblem, which is only a matter of determining x . Recall that we defined the polyp popu-lation density P in (3), and let b = ( b , b , . . . , b d ). Then the equation for x is given by x = λ ( b · x ) ϕ ( P ) . a = ( a , a , . . . , a d ). This immediately implies the identities x = x a , P = x Ω d (cid:88) k =2 p k a k , and b · x = ( b · a ) x . Altogether, this shows that a vector x = ( x , . . . , x d ) is a fixed point for the map f ( λ, · ) ifand only if x = x a and its first component x satisfies the nonlinear equation x = λ ( b · a ) x ϕ (cid:32) x Ω d (cid:88) k =2 p k a k (cid:33) . (7)From this equation, one can then determine all fixed points of the coral population model.Notice that we clearly have the trivial solution x = 0 for all values of the parameter λ , whichcorresponds to an extinct population. An important biological parameter for the coral population is the total number of larvaeproduced by a single colony in its entire life span. This number only depends on the birthand survival rates, and one can easily see that it is given by R = λb + λb S + λb S S + · · · + λb d S d − S d − . . . S = λ d (cid:88) i =1 a i b i . (8)The number R is called the basic reproduction number . Using the notation from the lastsubsection, the above equation can be rewritten as R = ( b · a ) λ . (9)In particular, while it is possible to vary R in such a way that the relationship between thebirth rate constants vary, under our assumptions, the vectors b and a are fixed constantvectors, and we therefore have a fixed linear relationship between R and λ . To make it easyto compare our results with those of previous papers, we have chosen to plot all bifurcationdiagrams with respect to the basic reproduction number R . We now turn our attention to a description of the bifurcation diagram of the fixed points forthe coral population system. This diagram is shown in Figure 4, where the set of fixed pointsis plotted in terms of the reproductive number R versus polyp population density P . Thecolor in the diagram depicts the stability of the fixed points, and the diagram indicates theexistence of three bifurcation points: a saddle-node and a Neimark-Sacker bifurcation on thenontrivial branch, which itself bifurcates from the trivial branch at a transcritical bifurcation.While subsequent sections of this paper will be used to verify the bifurcation diagram usingcomputer-assisted proofs, the remainder of the current subsection is devoted to the discussionof dynamical aspects which are observed through numerical simulations.7igure 4: The bifurcation diagram of polyp density P as a function of the reproductivenumber R . While the diagram covers the range R ∈ (12 , R ≈ d = 13 age groups. The bifurcationdiagram in Figure 4 was computed using a numerical continuation method starting at repro-duction number R = 300, and allowing R to decrease. There appears to be a saddle-nodepoint for R ≈ .
28 (which corresponds to λ ≈ . R of the fixed points begins to increase again. In Section 4 we use a computer-assistedproof to rigorously validate this saddle-node bifurcation point. The curve continues furtheruntil the population density reaches zero, which corresponds to an extinct population. Wewill see later that the extinction point can be found explicitly, and that it occurs at R ≈ . λ ≈ . x = 0 canreadily be determined from the Jacobian matrix of f at the origin, and this shows that theextinction fixed point is stable for small R , corresponding to low fitness, and unstable for alllarger values of the basic reproduction number R , with instability index 1. The bifurcationbetween the extinction fixed point being stable and unstable occurs at the transcritical bi-furcation point. All of these statements will be established rigorously in Section 4, includingthe appearance of the transcritical bifurcation point. Unlike the other two bifurcation points,no computer-assisted proofs are necessary along the trivial solution.As mentioned before, the stability of the fixed points x ∗ ∈ R is indicated by color, withblue indicating stable fixed points and red representing unstable ones. The local stability ateach fixed point in Figure 4 is determined numerically, based on whether all the eigenvaluesof the Jacobian matrix D x f ( λ, x ) lie inside the unit circle or not. In the bifurcation diagram,we have not distinguished the index of the stability. If at least one of the eigenvalues liesoutside the complex unit circle, then the fixed point is colored red, meaning unstable.8 .5 Oscillations Figure 4 only shows the existence and stability behavior of fixed point solutions. But whatabout the dynamical behavior of the system? In this last subsection of Section 2, we focus ondynamical aspects of the model, in particular its oscillatory behavior on attracting invariantcircles that form as a result of the Neimark-Sacker bifurcation. For a fixed parameter value
R > . y ∈ R . Atreproduction number R = 8 .
744 (which corresponds to λ = 0 . R = 29 .
15 (correspondingto λ = 1), if we start at initial conditions ranging roughly from 0 . y to 2 y , where y is avector of age-structured initial number of colonies which was chosen with polyp populationdensity P = 1500, then solutions converge to a nontrivial stable fixed point. There is also anunstable fixed point denoted by the red line. In addition, one can observe bistability at thisparameter value. If we start at a smaller value of P , such as for example at initial populationswith polyp population density smaller than 0 . y , solutions converge to zero, i.e., the coralpopulation becomes extinct. At the basic reproduction number R = 87 . λ = 3), thoughit takes longer time than 100 years, the solutions still converge to a stable nontrivial fixedpoint. In contrast, at R = 160 .
31 ( λ = 5 . P = 1 . y oscillate. Weused connected lines to show these oscillations more effectively, but recall that the map is infact discrete.The oscillations seen in the lower right subplot of Figure 5 form as a result of the Neimark-Sacker bifurcation. The fixed point stability switches from stable to unstable, and an invariantcircle gains stability. Trajectories with initial conditions near fixed points but after thebifurcation are displayed in Figure 6. Perturbations around an unstable fixed point arerepelled from the fixed point after the bifurcation, converging to an invariant closed curve.As the parameters R and λ increase, the size of the closed curve also increases, and theminimum population of a curve approaches the extinction point at the origin. That is,red coral populations become vulnerable at a large reproduction number, and a very smallperturbation of the population would endanger the survival of the population despite theexisting long recovery cycle.In order to better understand the stable invariant limit cycles that form after bifurcation,we have computed the rotation number, meaning the average angle of rotation per iterate, as afunction of the parameter R . Specifically, we used the projection to the x x -plane to computethe rotation numbers. Our computations are performed using the weighted Birkhoff averagemethod described in [7]. Figure 7 shows cycles at a ten distinct parameter values on the left,and for 500 distinct parameters on the right. The corresponding rotation numbers are shownin Figure 8. The values are angles, but they are rescaled to have values in the range (0 , , P , and they aresimulated over a time frame of 100 years each, at various parameter values.for a series of test parameters. In these test parameters, the answer differs by 10 − or less.Note that we would expect to see a devil’s staircase in the rotation numbers at theparameter values when there are periodic orbits, but what we see looks smooth even whenquite zoomed in. This is due to the fact that the periodic orbits are extremely high period.In particular, we are able to use a Farey tree calculation to find the smallest denominator,corresponding to the lowest period, of a periodic orbit for the case of a rational rotationnumber for this range of rotation numbers, using the method in [2, 14]. In particular, we findthat the lowest denominator in the range [0 . , . / R . That is, for each point in theinvariant circle, we graph how much the population is changing in one iterate (corresponding10igure 6: After the Neimark-Sacker bifurcation, oscillating orbits appear. After removingtransients in the orbit, the orbit lies on an invariant closed curve. On the left, we plot the x -and x -components of these limit cycles. As the parameter R increases, the size of the closedcurve increases. For large values of R , the coral population is close to the extinction pointat the origin. On the right, the same orbits are shown with respect to R , along with thecorresponding unstable fixed points at the same parameter value.Figure 7: Top: Invariant cycles for ten (left) and 500 (right) different parameter values. Eventhough we are guaranteed that some of the cycles contain stable periodic orbits, the periodsare sufficiently high and the parameter ranges for which they exist are sufficiently small thatit is hard to see them even in a close zoom (not depicted). Each orbit was computed using100,000 iterates.to one year) at each point in the invariant circle. The smallest angle difference, correspondingto the slowest change, occurs for angle ≈ . x , x ), computed with respect to the point( x , x ) = (2500 , R values. The minimum occurs at the angle pointing towards theextinction point.angle differences are very close to zero. Thus the population remains extremely vulnerablefor a particularly long time. We now turn to the rigorous validation of fixed points, both for regular and bifurcationvalues. Our general approach is the constructive implicit function theorem from [15]. Thisis a rigorous result that combines with a numerical interval arithmetic calculation to giverise to a validated method for finding a branch in the zero set of a function which dependson a single parameter. In the following four subsections, we will first recall the constructiveimplicit function theorem, and then define an extended system which can be used for pseudo-arclength continuation. After that, we prove two results which form the basis of our approach,and describe the necessary preconditioning for the coral population model application.
Before stating the full result, here is a summary. Given an approximate zero ( α ∗ , x ∗ ) of afunction G ( α, x ) where x is contained in a Banach space and α ∈ R , under certain hypotheseson G and its derivatives evaluated at the approximate zero ( α ∗ , x ∗ ), combined with Lipschitzestimates near this point, there exist two regions in parameter and phase space. First, the accuracy region , which contains a curve of the zero set. Second, a uniqueness region , in whichthat zero set curve is unique. See the schematic in Figure 9. The blue dot shows the initialapproximate zero. The orange curve is the zero set curve, which is guaranteed to lie withinthe accuracy region (the blue region). Note that the approximate zero does not in generallie on the zero set. The accuracy region is contained within the uniqueness region, shown inorange. The uniqueness region is largest in phase space when the parameter is closest α ∗ . As12 u (cid:1) α UniquenessAccuracy
Figure 9: A schematic depiction of the constructive implicit function theorem. The theoremguarantees that under appropriate hypothesis, an approximate zero (blue dot) guaranteesthat within a uniqueness region (orange region) there is a curve in the zero set with a uniquepoint at each fixed α value (red curve), and the this curve is located within an accuracyregion (blue region). The uniqueness region contains the accuracy region. It is bounded innorm by straight lines, and the accuracy region is bounded in norm by parabolas.the parameter varies, the uniqueness region shrinks (meaning we have worse isolation). Theconstructive implicit function theorem guarantees that the uniqueness region is characterizedby a linear norm condition, as depicted by the straight sides in the schematic diagram. Theaccuracy region has best (i.e., smallest) accuracy when the parameter is near the parameterof the original point α ∗ . The accuracy region grows (meaning we have worse accuracy) witha quadratic norm condition. This is depicted schematically by its parabolic shape. We nowstate the formal theorem. Theorem 3.1 (Constructive Implicit Function Theorem) . Let P , X , and Y be Banach spaces,suppose that the nonlinear operator G : P × X → Y is Fr´echet differentiable, and assume thefollowing hypotheses.(H1) Small residual: There exists a pair ( α ∗ , x ∗ ) ∈ P × X and a (cid:37) > such that (cid:107) G ( α ∗ , x ∗ ) (cid:107) Y ≤ (cid:37) . (H2) Bounded derivative inverse: There exists a constant K > such that (cid:13)(cid:13) D x G ( α ∗ , x ∗ ) − (cid:13)(cid:13) L ( Y , X ) ≤ K , where (cid:107) · (cid:107) L ( Y , X ) denotes the operator norm in L ( Y , X ) .(H3) Lipschitz bound: There exist positive real constants L , L , (cid:96) x , and (cid:96) α ≥ such thatfor all pairs ( α, x ) ∈ P × X with (cid:107) x − x ∗ (cid:107) X ≤ (cid:96) x and (cid:107) α − α ∗ (cid:107) P ≤ (cid:96) α we have (cid:107) D x G ( α, x ) − D x G ( α ∗ , x ∗ ) (cid:107) L ( X , Y ) ≤ L (cid:107) x − x ∗ (cid:107) X + L (cid:107) α − α ∗ (cid:107) P . H4) Lipschitz-type bound: There exist positive real constants L and L , such that for allparameters α ∈ P with (cid:107) α − α ∗ (cid:107) P ≤ (cid:96) α one has (cid:107) D α G ( α, x ∗ ) (cid:107) L ( P , Y ) ≤ L + L (cid:107) α − α ∗ (cid:107) P , where (cid:96) α is the constant that was chosen in (H3).Finally, suppose that K (cid:37)L < and K(cid:37) < (cid:96) x . (10) Then there exist pairs of constants ( δ α , δ x ) with ≤ δ α ≤ (cid:96) α and < δ x ≤ (cid:96) x , as well as KL δ x + 2 KL δ α ≤ and K(cid:37) + 2 KL δ α + 2 KL δ α ≤ δ x , (11) and for each such pair the following holds. For every α ∈ P with (cid:107) α − α ∗ (cid:107) P ≤ δ α there existsa uniquely determined element x ( α ) ∈ X with (cid:107) x ( α ) − x ∗ (cid:107) X ≤ δ x such that G ( α, x ( α )) = 0 .In other words, if we define B X δ = { ξ ∈ X : (cid:107) ξ − x ∗ (cid:107) X ≤ δ } and B P δ = { p ∈ P : (cid:107) p − α ∗ (cid:107) P ≤ δ } , then all points of the solution set of the equation G ( α, x ) = 0 in the set B P δ α × B X δ x lie on thegraph of the function α (cid:55)→ x ( α ) . In its classical form, the implicit function theorem is one of the central tools of bifurcationtheory. Not only can it be used to establish the existence of small solution branches innonlinear parameter-dependent equations, but by applying it as a tool to modified problemsit can frequently be used to provide sufficient conditions for bifurcations. For example, thecelebrated Crandall-Rabinowitz result [6] on bifurcation from a simple eigenvalue proves theexistence of a bifurcating branch by applying the implicit function theorem to a modificationof the original nonlinear problem which removes the trivial solution. The constructive implicitfunction theorem can similarly be used as a tool for bifurcation analysis, yet in a computer-assisted proof setting. In fact, some first applications in this direction have already beenprovided in [12, 15]. With the current paper, we add two more applications.More precisely, in the following we will be applying Theorem 3.1 in two different situations.In the remainder of this section, we apply it for branches of regular points. Through theintroduction of a suitable extended system we can reformulate a validated step of pseudo-arclength continuation as an application of the constructive implicit function theorem to thisextended system. Combined with suitable linking conditions, this establishes the existenceof entire branches covered by slanted boxes.In addition, in Section 4 we use Theorem 3.1 to validate bifurcation points. In thatsetting, and motivated by our earlier work [12], we will apply the theorem to an extendedsystem without any parameter, as the parameter will be incorporated into the function forwhich we find a root. This parameter-free case means that we no longer need to find theLipschitz constants relevant to the parameter variations, and we set these unused constantsequal to zero. 14
50 100 150 200 250 300 R P R un i qu e n e ss Figure 10: Left, the validated bifurcation diagram of polyp density P as a function of thereproductive number R , along with the three validated bifurcation points. The blue curveconsists of 5000 continuation steps, corresponding to 5000 linked boxes, for the precondi-tioned map with α = 0 . δ α . The initial validated box contains ( R, P ) = (300 , R, P ) = (71 . , . To elaborate further on the validation of regular fixed points, the constructive implicit func-tion theorem as stated in [15] only applies to a single region, validated at a single point. Thesame paper contains a version of this theorem for slanted boxes, using natural continuationin order to validate a branch of solutions by linking their validation sets to validate a largerportion of the branch. However, natural continuation leaves something to be desired in termsof efficiency. In this section, we develop a method of validation of bifurcation branches us-ing pseudo-arclength continuation which allows for the direct application of the constructiveimplicit function theorem, and apart from Lipschitz estimates, only requires estimates at asingle point in each box. This method is an improvement on the previous natural continua-tion method in that we can continue at limit points without having to change coordinates.The methods in this section apply for regular orbits along branches. In the next section,we will show how to adapt the constructive implicit function theorem in order to rigorouslyvalidate bifurcation points.Before launching into further technicalities, we describe our results. Applying the pseudo-arclength continuation method to a preconditioned version of the coral model (preconditioningis discussed in Section 3.4 below), the resulting rigorously validated curve of fixed points isshown in Figure 10. While Figure 4 shows a similar picture, the distinction is that thosepoints were found using numerical methods, and though we have a priori error estimates for15hese methods, we cannot guarantee existence or accuracy. In contrast, the points shownon the new figure are rigorously validated. The depicted points are an accurate indicationof existing fixed points of the system, with known and validated accuracy and uniquenessregion. In particular, the accuracy of our solutions is known individually for each separatebox, and is always less than 1 . · − , where the error in x ∈ R is measured in themaximum norm. Figure 10 shows the norm of the uniqueness for each separate box. Theuniqueness shrinks when the curve approaches zero. This is not surprising, since x = 0 ispart of the zero set, putting a barrier on the size of the uniqueness region.We now proceed with the constructive implicit function theorem for a validated pseudo-arclength continuation. In each continuation step we use continuation in a box with slantedsides, where the predictor step is performed along the middle of the box in the direction aspecified vector ( µ, v ) (usually the estimated tangent to the zero set curve), and the correctorstep uses a computation such as Newton’s method to refine the estimate. This refinementis performed in a direction orthogonal to the predictor direction ( µ, v ). This is depicted inFigure 11. The left-hand image is a schematic diagram showing the box with its midlinebetween two blue dots. The midline is the estimated tangent line in the direction ( µ, v ). Ourvalidation gives us a maximum length of the box for which we can guarantee accuracy anduniqueness of the solution. The predictor, shown with a red dot, must be chosen inside thatbox. The corrector, shown with a green dot is along an orthogonal line to the midline. Theright-hand image shows the accuracy region in blue and the uniqueness region in orange.Note that the uniqueness region has large width near the starting point, and the accuracyregion grows towards the ending point. In Figure 11, the uniqueness region for the box isapproximately diamond shaped, whereas in Figure 12, the box is not only slanted but alsohas a uniqueness region which is asymmetric, more of a half-diamond. The half-diamondshape is in fact only half of the uniqueness box. In particular, as we are merely continuing inone direction, which in Figure 12 is to the left, we only show one side of the uniqueness box.The fact that we could continue to the right as well is not relevant for our continuation.We now turn to the technical details of this approach. For this, let F : R × U → U ,where U denotes an arbitrary Euclidean space. Our goal is to implement pseudo-arclengthcontinuation based on Theorem 3.1 to find branches of zeros of the nonlinear function F . Forthe specific application of this paper, we will consider U = R and F ( λ, x ) = f ( λ, x ) − x ,where f is the coral model. Nevertheless, we use the more general notation based on F to indicate that these methods are general. In fact, the methods readily generalize to theBanach space setting as well. However, in this paper for convenience of notation we onlyconsider the Euclidean space case. For any ( λ , u ) ∈ R × U , an approximate zero of F , andfor a fixed direction vector ( µ , v ) ∈ R × U , define G : R × ( R × U ) → R × U as follows G ( α, ( σ, x )) = (cid:32) µ σ + v t xF ( λ + αµ + σ, u + αv + x ) (cid:33) . (12)The zeros of G as the parameter α varies correspond to the pseudo-arclength continuationsolutions of F for a single continuation box. The first component of the function G guaranteesthat the pair ( σ, x ) is orthogonal to the direction ( µ , v ). As we will show in the nextsubsection, one can apply the constructive implicit function theorem from [15] directly to theextended function G and thereby perform rigorously validated pseudo-arclength continuation.16 λ ∗ k , u ∗ k )( λ ∗ k , u ∗ k ) + α ( µ k , v k ) + ( σ, x )( λ ∗ k , u ∗ k ) + δ α ( µ k , v k ) ( σ, x ) ( λ ∗ k +1 , u ∗ k +1 )( λ ∗ k , u ∗ k )( λ ∗ k , u ∗ k ) + α ∗ ( µ k , v k ) ( λ ∗ k , u ∗ k )( λ ∗ k , u ∗ k ) + δ α ( µ k , v k ) Figure 11: A schematic diagram of the the pseudo-arclength continuation method. Leftimage: The result guarantees a uniqueness region for the zero set. This takes place in anadapted coordinate system, meaning that the box is slanted, but the uniqueness region is stillbounded by straight lines. Since we only continue the curve in one direction, this figure onlydepicts the left half of the uniqueness region. The center line segment of this region is givenby ( λ ∗ k , u ∗ k ) + α ( µ k , v k ) for 0 ≤ α ≤ δ α . At a fixed α value, we use Newton’s method to findthe next approximate zero along the line ( λ ∗ k , u ∗ k ) + α ( µ k , v k ) + ( σ, x ), where ( σ, x ) denotesthe vector pointing from ( λ ∗ k , u ∗ k ) + α ( µ k , v k ) to the point ( λ ∗ k , u ∗ k ) + α ( µ k , v k ) + ( σ, x ), andwhich is orthogonal to ( µ k , v k ). Middle image: After we fixed the value α = α ∗ , we labelthis next approximation ( λ ∗ k +1 , u ∗ k +1 ). Right image: Inside the uniqueness region (orange) isan accuracy region (blue). The accuracy region is bounded by curves which are parabolic innorm in the adapted coordinate system.Since we will need them later, we close this subsection by explicitly stating the derivativesof G with respect to both the variables ( σ, x ) and with respect to the parameter α . Theseare respectively given by D ( σ,x ) G ( α, ( σ, x )) = (13) (cid:32) µ v t D λ F ( λ + αµ + σ, u + αv + x ) D u F ( λ + αµ + σ, u + αv + x ) (cid:33) , as well as D α G ( α, ( σ, x )) = (14) (cid:32) D λ F ( λ + αµ + σ, u + αv + x ) µ + D u F ( λ + αµ + σ, u + αv + x ) v (cid:33) . We are now in a position to start establishing assumptions under which we can validatea branch in the zero set of F using pseudo-arclength continuation. For this we need thefollowing modified set of assumptions. For the purposes of this paper, we use the vectornorm (cid:107) ( α, x ) (cid:107) = max {| α | , (cid:107) x (cid:107) U } for all ( α, x ) ∈ R × U , even though this could easily bemodified.(P1) We assume both (cid:107) F ( λ , u ) (cid:107) U ≤ (cid:37) and (cid:107) D λ F ( λ , u ) µ + D u F ( λ , u ) v (cid:107) U ≤ ξ . (15)17P2) Assume that there exists an explicit constant K > D ( σ,x ) G (0 , (0 , (cid:32) µ v t D λ F ( λ , u ) D u F ( λ , u ) (cid:33) , i.e., we suppose that (cid:13)(cid:13) D ( σ,x ) G (0 , (0 , − (cid:13)(cid:13) L ( R × U, R × U ) ≤ K .
For this, we interpret the matrix as a linear map on the product space R × U , and theoperator norm is the norm in L ( R × U, R × U ).(P3) Let M , M , M , and M be Lipschitz constants such that for all pairs ( λ, u ) whichsatisfy (cid:107) u − u (cid:107) ≤ d u and | λ − λ | ≤ d λ we have the estimates (cid:107) D u F ( λ, u ) − D u F ( λ , u ) (cid:107) L ( U,U ) ≤ M (cid:107) u − u (cid:107) U + M | λ − λ | , (cid:107) D λ F ( λ, u ) − D λ F ( λ , u ) (cid:107) L ( R ,U ) ≤ M (cid:107) u − u (cid:107) U + M | λ − λ | , where as usual we will identify the norm in L ( R , U ) with the norm (cid:107)·(cid:107) U in the following.We would like to point out that all of the above three conditions are formulated in terms ofthe nonlinear parameter-dependent function F and an approximate solution ( λ , u ) of theequation F ( λ, u ) = 0.We now turn our attention to the extended system described by the operator G intro-duced in (12). It turns out that the above three assumptions are tailor-made to establishthe hypotheses (H1) through (H4) from the constructive implicit function theorem for themapping G . One can easily see that (P1) implies (cid:107) G (0 , (0 , (cid:107) R × U ≤ (cid:37) , i.e., hypothesis (H1) is satisfied. Furthermore, using the explicit derivative formulas from theend of the last subsection, the assumption (P2) immediately yields the estimate (cid:107) D ( σ,x ) G (0 , (0 , (cid:107) L ( R × U, R × U ) ≤ K , which establishes (H2). It remains to show that (P3) furnishes the estimates in (H3) and (H4).For this, let ξ be defined as in (15), and define the four constants L = max( M + M , M + M ) ,L = ( M + M ) (cid:107) v (cid:107) U + ( M + M ) | µ | ,L = ξ ,L = ( M (cid:107) v (cid:107) U + M | µ | ) (cid:107) v (cid:107) U + ( M (cid:107) v (cid:107) U + M | µ | ) | µ | . L through L are the Lipschitz constants for the extended function G as required by (H3) and (H4). For this, first note that in view of (13) we have D ( σ,x ) G ( α, ( σ, x )) − D ( σ,x ) G (0 , (0 , (cid:32) D λ F ( w ) − D λ F ( w ) D u F ( w ) − D u F ( w ) (cid:33) , where D λ F and D u F are evaluated at w = ( λ + αµ + σ, u + αv + x ) and w = ( λ , u ).Then one can readily see that (H3) follows from (P3) and the estimates (cid:107) D ( σ,x ) G ( α, ( σ, x )) − D ( σ,x ) G (0 , (0 , (cid:107) L ( R × U, R × U ) ≤ (cid:107) D u F ( λ + αµ + σ, u + αv + x ) − D u F ( λ , u ) (cid:107) L ( U,U ) + (cid:107) D λ F ( λ + αµ + σ, u + αv + x ) − D λ F ( λ , u ) (cid:107) L ( R ,U ) ≤ M ( | α |(cid:107) v (cid:107) U + (cid:107) x (cid:107) U ) + M ( | α || µ | + | σ | )+ M ( | α |(cid:107) v (cid:107) U + (cid:107) x (cid:107) U ) + M ( | α || µ | + | σ | )= ( M + M ) (cid:107) x (cid:107) U + ( M + M ) | σ | + (( M + M ) (cid:107) v (cid:107) U + ( M + M ) | µ | ) | α | = L (cid:107) ( σ, x ) (cid:107) R × U + L | α | . Similarly, using (14) one can show that (H4) follows from (P1) and (P3), in combination withthe inequalities (cid:107) D α G ( α, (0 , (cid:107) L ( R , R × U ) ≤ (cid:107) D λ F ( λ , u ) µ + D u F ( λ , u ) v (cid:107) U + (cid:107) D u F ( λ + αµ , u + αv ) v − D u F ( λ , u ) v (cid:107) U + (cid:107) D λ F ( λ + αµ , u + αv ) µ − D λ F ( λ , u ) µ (cid:107) U ≤ ξ + ( M (cid:107) v (cid:107) U + M | µ | ) | α |(cid:107) v (cid:107) U + ( M (cid:107) v (cid:107) U + M | µ | ) | α || µ | = L + L | α | . Altogether, these estimates lead to the following result.
Theorem 3.2 (Pseudo-arclength continuation for a branch segment) . Consider the fixedpairs ( u , λ ) and ( v , µ ) in R × U , let d λ and d u be two positive constants, and suppose thatour hypotheses (P1), (P2), and (P3) are satisfied. Moreover, assume that both K (cid:37) < and K(cid:37) < d u hold. Then we can choose constants < δ α ≤ d λ , < δ u ≤ d u , where δ α (cid:107) ( µ , v ) (cid:107) + δ u ≤ min( d u , d λ ) , and such that KL δ u + 2 KL δ α ≤ and K(cid:37) + 2 KL δ α + 2 KL δ α ≤ δ u . λ ∗ k +1 , u ∗ k +1 )( λ ∗ k , u ∗ k ) -- ( λ ∗ k , u ∗ k )( λ ∗ k +1 , u ∗ k +1 )( λ ∗ k +1 , u ∗ k +1 ) + ( σ new , x new ) Figure 12: Left image: Associated with each successive approximation, there is a uniquenessregion and an accuracy region. Right image: In order to guarantee that the k -th and ( k + 1)-st region enclose the same component of the zero set (the green curve), we must verify thelinking condition. This requires that the accuracy curve of the ( k + 1)-st box at α = 0 (suchas the blue point on the upper edge of the ( k + 1)-st blue box) is contained in the uniquenessregion of the k -th box (orange region). Then for every α ≤ δ α there exists a unique ( σ, x ) in the zero set of G with (cid:107) ( σ, x ) (cid:107) ≤ δ u .These statements guarantee that there is a unique element of the zero set of F whichlies on the hyperplane orthogonal to the center line in the slanted box between ( λ , u ) and ( λ + δ α µ , u + δ α v ) and passes through the point ( λ + αµ , u + αv ) . This unique zero isgiven by ( λ + αµ , u + αv ) + ( σ, x ) . Additionally, let δ min = 2 K(cid:37) .
Then for α = 0 we can guarantee that the resulting pair in the zero of G is accurate within δ min of ( λ , u ) , and this zero is unique within the set (cid:107) ( σ, x ) (cid:107) ≤ min { (2 KL ) − , d u , d α } .Proof. To show the theorem we follow the proof of [15, Theorem 5]. Aside from the changesin the Lipschitz constants which have already been derived before the formulation of thetheorem, the only changes to the cited proof are due to the fact that for a fixed parameterof G , the values of both the parameter λ and the phase space value x of F can vary. Therefore,in order to guarantee that the Lipschitz estimates on F hold, we need to assure that for every α ≤ δ α and all (cid:107) ( σ, x ) (cid:107) ≤ δ u the norm (cid:107) α ( µ , v ) + ( σ, x ) (cid:107) is bounded by both d u and d λ .This immediately leads to the additional constraints in the formulation of the theorem.The above theorem gives a method for validating a branch segment of the zero set withina single slanted box. In practice we use this result successively to validate a whole solutionbranch. For each pair ( λ ∗ k , u ∗ k ), and for the approximate tangent ( µ k , v k ), we then define anextended function G k , and validate a branch segment for F within the k-th box. For a fixedparameter value α k ≤ δ α , we then use Newton’s method to find an approximate zero of F which is orthogonal to ( µ k , v k ), i.e., which is a zero of G k . We abbreviate this approximatezero as ( λ ∗ k +1 , u ∗ k +1 ), and can now repeat the entire process for the ( k + 1)-st branch segment,see also Figure 12. What remains to be shown is that the successive validated boxes arelinked, meaning that the branch segment in the k -th box and the branch segment in the( k + 1)-st box are on the same branch. That is, the accuracy region of the ( k + 1)-st box has20o be contained within the uniqueness region of the k -th box at the point α k where we madethe numerical estimate. We give the linking condition for two boxes in the next theorem. Theorem 3.3 (Linking branch segments) . Let δ k +1 , min = 2 K k +1 (cid:37) k +1 be the accuracy of thesolution ( λ ∗ k +1 , u ∗ k +1 ) = ( λ ∗ k + α k µ k + σ ∗ , u ∗ k + α k v k + x ∗ ) . In order to guarantee that the two validated boxes are linked, we require the estimates | α k | + δ k +1 , min (cid:107) ( µ k , v k ) (cid:107) < δ k,α and | ( σ ∗ , x ∗ ) | + δ k +1 , min < δ k,u . Proof.
The accuracy of the ( k + 1)-st solution at α = 0 is given by δ k +1 , min . That is, thereexists a unique exact solution to F = 0 of the form(˜ λ, ˜ u ) = ( λ ∗ k +1 + σ new , u ∗ k +1 + x new ) , where (cid:107) ( σ new , x new ) (cid:107) < δ k +1 , min . In order to derive our linking condition we need to establishthat this solution is contained in the uniqueness region of the k -th segment. We can thereforewrite (˜ λ, ˜ u ) − ( λ ∗ k , u ∗ k ) = ( α k + α + )( µ k , v k ) + ( σ ∗ + σ + , x ∗ + x + ) , where ( σ new , x new ) = α + ( µ k , v k ) + ( σ + , x + ), and ( µ k , v k ) is orthogonal to ( σ + , x + ). Thus wehave (cid:107) α + ( µ k , v k ) + ( σ + , x + ) (cid:107) < δ k +1 , min . By the orthogonality of the two vectors, both the estimate | α + |(cid:107) ( µ k , v k ) (cid:107) < δ k +1 , min and theestimate (cid:107) ( σ + , x + ) (cid:107) < δ k +1 , min are satisfied. In order to satisfy the linking condition, we haveto require that both | α k + α + | < δ k,α and (cid:107) ( σ ∗ + σ + , x ∗ + x + ) (cid:107) < δ k,u hold. This translatesinto the conditions | α k + α + | ≤ | α k | + δ k +1 , min (cid:107) ( µ k , v k ) (cid:107) < δ k,α , as well as (cid:107) ( σ ∗ + σ + , x ∗ + x + ) (cid:107) ≤ (cid:107) ( σ ∗ , x ∗ ) (cid:107) + δ k +1 , min < δ k,u . This completes the proof of the theorem.
If we use the above method on the coral system, it is extremely slow to produce the bifurcationdiagram. This is due to the different relative sizes of the components of the population andthe parameter. We are able to significantly speed up the method by using preconditioning.In particular, for k = 1 , . . . , d let˜ f k ( ˜ R, ˜ u ) = f k (100 ˜ R, ( s ˜ u , . . . , s d ˜ u d )) s k , where s , . . . , s d are empirically determined positive scale constants. Then it is clear that ifwe write ( R, u ) = (100 ˜ R, ( s ˜ u , . . . , s d ˜ u d )), then ( R, u ) is a fixed point of f if and only if ( ˜ R, ˜ u )21s a fixed point of the preconditioned map ˜ f . However, the map ˜ f is better scaled in thesense that we expect all components and the parameter to be of the same order of magnitude.Therefore the pseudo-arclength continuation can be performed more efficiently. In particular,we find that the size of δ α in the preconditioned version is (in comparable coordinates) aroundtwo orders of magnitude larger than those for the system without modification. This meansthat we are able to validate a much larger portion of the bifurcation diagram with the samenumber of continuation steps. Figure 10 shows 5000 continuation steps for the preconditionedcase starting at R = 300 in the upper right corner, shown in blue. For comparison purposes,4000 continuation steps are shown in red for the unmodified case. The bifurcation curve goesthrough a limit point and almost to (cid:107) u (cid:107) = 0 for the preconditioned case, but is hardly evena visible piece of red curve for the original unmodified map. A similar preconditioning isperformed in the case of the bifurcation points, as described in the next section. In this section, we discuss the validation of the bifurcation points. Namely, we have used acomputer-assisted proof to validate the Neimark-Sacker bifurcation point, where the invariantcircles form in Section 4.1 and the saddle-node bifurcation point in Section 4.2. In eachcase, to do so we create an extended system H such that H = 0 guarantees the neededconditions for a bifurcation point. We then apply the constructive implicit function theoremto H . In both cases, we use interval arithmetic for a separate computational validationof the extra transversality and nondegeneracy conditions. We also prove that there is atranscritical bifurcation point on the extinction axis. However, this last case does not requirea computer-assisted proof for validation, since the calculations are simple enough for a closedform calculation. In Sections 2.4 and 2.5, we observed that at (
R, P ) ≈ (154 . , R > .
1, typical initial conditions converge to popu-lations which are oscillating in time. This is the behavior associated with a Neimark-Sackerbifurcation. In this section we detail the process of rigorous validation of the Neimark-Sackerbifurcation point seen in the upper right corner of Figure 4. While this is the first timethat a rigorous validation of a Neimark-Sacker bifurcation has been performed in this way,rigorous validation of Hopf bifurcations was performed in [18] in the context of ordinary andpartial differential equations, but using a quite different method. Rather than consideringconditions along a curve of fixed points or equilibria, instead the method used a validatedcontinuation of periodic orbits with a renormalization technique, validating that there wasa bifurcation of equilibria at the turning point of this invariant closed curve of solutions.Moreover, computer-assisted proofs were used in [5] to rigorously establish an invariant circlein a two-dimensional map, which is created via a Neimark-Sacker bifurcation. They do not,however, establish the bifurcation point itself directly. While it would be interesting to adapttheir method to the coral model, this lies beyond the scope of the current paper.22e now proceed with our validation of the Niemark-Sacker bifurcation. As a first step,we state the standard theoretical Neimark-Sacker bifurcation theorem found in a bifurca-tion theory textbook. We then show how to adapt this classical result to create a rigorouscomputer-assisted bifurcation theorem.
Theorem 4.1 (Neimark-Sacker bifurcation point) . There is a Neimark-Sacker bifurcationfor the coral system in (4) and (5) for the basic reproduction number R ∗ ≈ . and withpolyp population density P ∗ ≈ . The precise error bounds are stated in Table 2.
The remainder of this subsection is devoted to the proof of this theorem. Our approach isto verify the classical conditions for a
Neimark-Sacker bifurcation , as described for examplein [11] — and which we briefly review in the following. Consider a smooth map f : R × R d → R d . Furthermore, we begin by assuming the following two conditions:(a) Existence of a fixed point:
The map f has a fixed point at a specific parameter value,i.e., we assume that f ( λ , x ) = x .(b) Pair of imaginary eigenvalues on the unit circle:
The Jacobian matrix D x f ( λ , x )has exactly one simple conjugate pair of imaginary eigenvalues on the unit circle. Wedenote these eigenvalues by e ± iθ , for some angle 0 < θ < π .These two conditions have to be supplemented by another three transversality and nonde-generacy conditions, which will be stated in detail below. For this, however, we first need tointroduce some additional notation.Due to the implicit function theorem, as long as the Jacobian matrix in (b) does not havethe eigenvalue 1, there exists a smooth curve of locally unique fixed points, which we denoteby ( λ, x ( λ )). Moreover, we define A ( λ ) = D x f ( λ, x ( λ )) . We would like to point out that in our application to the coral system, the rigorously es-tablished existence of the branch of fixed points as a side effect also implies that along thebranch near the Neimark-Sacker point, the Jacobian matrix never has an eigenvalue 1.Now let p ∈ C d and q ∈ C d denote the right eigenvectors of A ( λ ) corresponding to e iθ and e − iθ , respectively, and normalized in such a way that (cid:104) p, q (cid:105) = 1, where the bracketnotation denotes the usual complex scalar product (cid:104) p, q (cid:105) := p t q . Finally, by Taylor’s formulawe can expand the function f in the form f ( λ , x ) − x = A ( λ ) x + 12 B ( x, x ) + 16 C ( x, x, x ) + O ( (cid:107) x (cid:107) ) , (16)where B and C denote the second- and third-order derivative terms at the point ( λ , x ) inthe form B i ( y, z ) = d (cid:88) j,k =1 ∂ f∂x j ∂x k ( λ , x ) y j z k and C i ( y, z, w ) = d (cid:88) j,k,l =1 ∂ f∂x j ∂x k ∂x l ( λ , x ) y j z k w l . After these preparations, we can now complete our description of the conditions needed forthe Neimark-Sacker theorem: 23 λ x P δ δ . .
286 1794 2689 1 . · − . · − (cid:37) K L (c) (d) (e)6 . · − .
000 4 . · . · − . − . · − Table 2: Validation constants for the system (17) at the Neimark-Sacker bifurcation point. Allvalues are written with four decimal places, unless less accuracy is known. For more efficientcomputation, we multiplied by a preconditioning matrix and determined the bounds (cid:37) , K ,and L . We selected a matrix close to the Jacobian matrix of H ns , whose inverse was used asa preconditioner. The accuracy constant δ and the isolation bound δ were derived using (cid:37) , K and L . For the three conditions (c), (d), and (e), which were checked separately afterthe validation involving H ns , we used an interval arithmetic enclosure of the approximatesolution with radius δ . Note that the angle in (d) is given in degrees.(c) Transversality condition:
Using the notation above, suppose thatRe (cid:18) e − iθ (cid:28) p, dAdλ ( λ ) q (cid:29)(cid:19) (cid:54) = 0 . (d) Nondegeneracy condition I:
Suppose that θ (cid:54) = π θ (cid:54) = 2 π . (e) Nondegeneracy condition II:
Suppose that Re ( e − iθ ( (cid:104) p, C ( q, q, ¯ q ) (cid:105) + 2 (cid:104) p, B ( q, ( I − A ) − B ( q, ¯ q )) (cid:105) + (cid:104) p, B ( q, ( e iθ I − A ) − B ( q, q )) (cid:105) ) ) (cid:54) = 0 . To summarize, the transversality condition implies that the pair of complex conjugate eigen-values at λ crosses the imaginary axis with nonzero speed. The first nondegeneracy conditionindicates that the eigenvalues e ± iθ are not k -th roots of unity for k = 1 , . . . ,
4. Since theproof of the Neimark-Sacker theorem is based on the Poincar´e normal form theorem, thiscondition excludes resonances. Finally, the left-hand side of the second nondegeneracy condi-tion gives the coefficient of the cubic term in the complex Poincar´e normal form, and its signdistinguishes between a sub- and super-critial Neimark-Sacker bifurcation. For more detailswe refer the reader to the part of [11, Section 5.4] devoted to the Neimark-Sacker bifurcation.Under the above conditions, the Neimark-Sacker theorem guarantees that a locally uniqueinvariant closed curve bifurcates from the set of fixed points at the point ( λ , x ). As alreadymentioned, the type of bifurcation depends on the sign of the left-hand side of (e).In order to create the validation version of this theorem, we use a suitable extended systemto validate assumptions (a) and (b). After having established an existence and uniquenessresult for this extended system, one can then validate conditions (c), (d), and (e) separatelyusing interval arithmetic. For convenience, we have converted the complex system into the24ollowing real system of equations. We are seeking zeros of the function H ns : R m → R m ,which is defined as H ns ( x, λ, w, u, a, b ) = f ( λ, x ) − xD x f ( λ, x ) w − aw + buD x f ( λ, x ) u − bw − aua + b − (cid:107) w (cid:107) − (cid:107) u (cid:107) − . (17)The first equation in the system is the fixed point condition. The second through fourthequations form the simple complex eigenvalue pair condition, where we write e ± iθ = a ± ib ,and the eigenvectors p and q are given by u ± iw , up to normalization. The last two equationsare included to single out a locally unique eigenvector.For a function of the form f : R × R d → R d , we have x ∈ R d , λ ∈ R , u, w ∈ R d , as well as a, b ∈ R . Therefore, the extended system H ns : R m → R m lives in dimension m = 3 d + 3. Inour numerical validation, we are working with a 13-dimensional system, implying that thisextended system has dimension 42.Using standard numerical methods, we obtained an approximate bifurcation point sat-isfying H ns ( x, λ, w, u, a, b ) = 0, for the function H ns in (17), and with values for R , λ , x ,and P as stated in Table 2. Since H ns is parameter free, we only seek rigorous solutionsof the extended system in (17) which satisfy H ns = 0 in R . Thus we only need to verifythe hypotheses of the constructive implicit function theorem which involve the values of (cid:37) , K , L , and (cid:96) x > (cid:37) and K by using interval arithmetic. While the bound (cid:37) canbe found in a straightforward way, the constant K cannot easily be found by using intervalarithmetic to compute matrix inverses. Therefore, we first compute an approximate numericalinverse. However, we still need a bound on the exact inverse, and a bound on the accuracy ofthe approximate inverse. This is required in both the computation of K and twice when weverify condition (e). The required quantities can be determined using the following lemma.While we apply this lemma only for matrices, it is stated for the case of Banach spaces. Lemma 4.2 (Inverse bounds) . Let A be a bounded linear operator between two Banach spaces,and let B be an approximate inverse of A . Assume further that (cid:107) I − BA (cid:107) ≤ (cid:37) < as well as (cid:107) B (cid:107) ≤ (cid:37) . Then A is one-to-one, onto, and we have both (cid:107) A − (cid:107) ≤ (cid:37) − (cid:37) and (cid:107) B − A − (cid:107) ≤ (cid:37) (cid:37) − (cid:37) . The bound on A − is due to a Neumann series argument, and the proof can be foundin [15]. In addition, the second bound is a consequence of (cid:107) B − A − (cid:107) ≤ (cid:107) I − BA (cid:107)(cid:107) A − (cid:107) .25aving described how the constants (cid:37) and K can be estimated rigorously, we now turnour attention to the Lipschitz constant L . It can be determined using the mean valuetheorem for multivariate functions from the calculations in (18) below. For this, supposethat the function H ns : R m → R m is differentiable and let h ij ( x ) = ( ∂ ( H ns ) i /∂x j )( x ). Then h ij : R m → R , and we let h : R m → R m × m denote the matrix-valued function with entries h ij .Throughout our computations, we used the maximum norms for vectors x , and the inducedmatrix norm for matrices A . Recall that one then has (cid:107) x (cid:107) = (cid:107) x (cid:107) ∞ = max i =1 ,...,m | x i | , as wellas (cid:107) A (cid:107) = (cid:107) A (cid:107) ∞ = max i =1 ,...,m (cid:80) mj =1 | A ij | . After these preparations, the mean value theoremimplies | h ij ( x ) − h ij ( y ) | ≤ max c ∈ D (cid:107)∇ h ij ( c ) (cid:107) (cid:107) x − y (cid:107) , where D denotes the line segment between the points x and y . Together with the definitionof the functions h ij one further obtains | h ij ( x ) − h ij ( y ) | ≤ max c ∈ D (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) ∂ ( H ns ) i ∂x ∂x j ( c ) , . . . , ∂ ( H ns ) i ∂x n ∂x j ( c ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) (cid:107) x − y (cid:107)≤ m max c ∈ D, k =1 ,...,m (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( H ns ) i ∂x k ∂x j ( c ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:107) x − y (cid:107) . This finally furnishes (cid:107) h ( x ) − h ( y ) (cid:107) = max i =1 ,...,m m (cid:88) j =1 | h ij ( x ) − h ij ( y ) |≤ max i =1 ,...,m m (cid:88) j =1 (cid:18) m max c ∈ D, k =1 ,...,m (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( H ns ) i ∂x k ∂x j ( c ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (cid:107) x − y (cid:107) . (18)The factor in front of (cid:107) x − y (cid:107) on the right-hand side is then the Lipschitz constant L , andit can be determined via interval arithmetic and automatic differentiation.Altogether, our rigorous computer-assisted proof of Theorem 4.1 can be summarized asfollows. After completing the validation of the conditions that guarantee that the constructiveimplicit function theorem holds, we are able to verify the accuracy and uniqueness regions forthe bifurcation point. In addition, we can use Intlab [13] to rigorously show that the Jacobianmatrix D x f ( λ , u ) has in fact only two eigenvalues on the unit circle, by verifying that theremaining eleven eigenvalues all lie inside the unit disk. This implies that a bifurcationoccurs within the specified error of the approximate bifurcation point. We then verify thatthis bifurcation is indeed a Neimark-Sacker bifurcation by showing that conditions (c), (d),and (e) hold using interval arithmetic on these conditions. Here are a few remarks which givea more detailed explanation: • For each condition, we show that the interval containing the exact answer does notcontain zero for (c) and (e), and does not contain any of the avoided angles for (d). • While we are able to work with real-valued quantities a, b, u, v in the initial calculationsof parts (a) and (b), we must switch to the complex case to verify the extra conditions26 λ x P δ δ .
28 0 . . . . · − . · − (cid:37) K L (c) (d)1 . · − . · − . − . · − Table 3: Validation constants for the extended system in (19) at the saddle-node bifurcationpoint. All values are written up to four decimal places. For more efficient computation, wemultiplied by a preconditioning matrix and obtained the bounds (cid:37) , K , and L . We selected amatrix close to the Jacobian matrix of H sn , whose inverse was used as a preconditioner. Theaccuracy constant δ and the isolation bound δ were derived using (cid:37) , K , and L . For thetwo conditions (c) and (d), which were checked separately after the validation involving H sn ,we used an interval arithmetic enclosure of the approximate solution with radius δ .(c), (d), (e), and we normalize the complex vectors p and q using the normalizationcondition (cid:104) p, q (cid:105) = 1. • We need to be able to guarantee that all three conditions are satisfied for the entireaccuracy region. Therefore we evaluate these conditions on an interval vector whosemidpoint is the approximate bifurcation point, and whose radius is δ . That is, everycomponent of the vector is an interval. The actual computed values of the conditions(c)-(e) are intervals, but the values given in Table 2 are the worst-case scenario values.Even with the interval calculations, conditions (c) and (d) are known to more than foursignificant digits, but condition (e) is only known to three digits of accuracy.This completes the proof of Theorem 4.1. In this section, we use a computer-assisted proof to show that there is a saddle-node bifur-cation point in the coral model. The precise result can be stated as follows.
Theorem 4.3 (Saddle-node bifurcation point) . The coral model in (4) and (5) has a saddle-node bifurcation point near the basic reproduction number R ∗ ≈ . , which corresponds tothe parameter value λ ∗ ≈ . , and for polyp population density P ∗ ≈ . . The preciseerror bounds are stated in Table 3. As in the previous subsection, the remainder of the present one is devoted to the veri-fication of this theorem via computer-assisted rigorous methods. In order to establish thetheorem, we need to verify the following conditions from the classical saddle-node bifurcationtheorem , see for example [11]. Let f : R × R d → R d be a smooth mapping. Furthermore,assume the following four conditions:(a) Existence of a fixed point:
The map f has a fixed point at a specific parameter value,i.e., we assume that f ( λ , x ) = x . 27b) Simple eigenvalue 1:
The Jacobian matrix D x f ( λ , x ) has a simple eigenvalue of 1.Let p and q denote the corresponding left and right eigenvectors, and suppose they arenormalized to satisfy p t q = 1.(c) Transversality condition:
Using the above notation we assume p t D λ f ( λ , x ) (cid:54) = 0 . (d) Nondegeneracy condition:
Now let A ( λ ) = D x f ( λ , x ), and consider the expansionof f given in (16). Then we suppose further that p t B ( q, q ) (cid:54) = 0 . Then the classical saddle-node bifurcation theorem guarantees a saddle-node bifurcation atthe pair ( λ , x ).In order to validate our bifurcation point using this theorem, we use again an extendedsystem of the form H sn = 0 to validate conditions (a) and (b), and then we verify condi-tions (c) and (d) separately afterwards. This time, the extended mapping H sn is a map H sn : R → R , and it is defined as H sn ( x, v, λ ) = f ( λ, x ) − xD x f ( λ, x ) v − v (cid:107) v (cid:107) − . (19)In order to validate (c), and (d), we use interval arithmetic for both of these conditions, andshow that 0 does not lie in the interval containing the resulting answer. Note that the vector q is just a multiple of v , and p can be found in a verified way using Intlab [13]. The summaryof the constants of this validation process is given in Table 3. This computer-assisted proofis quite similar to the one used for the Neimark-Sacker bifurcation in the last subsection,and therefore we do not give any more elaboration on the technique used to compute thesevalues. This completes the proof of Theorem 4.3. We close this section by showing that there is indeed a transcritical bifurcation on the trivialsolution curve, i.e., the extinction curve. This time, it is not necessary to perform a computer-assisted proof, as the bifurcation can be established directly by hand.
Theorem 4.4 (Transcritical bifurcation point) . For the coral population model in (4) and (5)there exists a transcritical bifurcation point for basic reproduction number R ∗ = c /c ≈ . ,which corresponds to the parameter value λ ∗ = R ∗ / ( b · a ) and to x ∗ = 0 ∈ R . Recall thatthe constants c and c were introduced in (2), and the vectors a and b were defined in (6)and the following paragraph. roof. It is clear from the model that x = 0 is a fixed point for all values of the parameter λ .Furthermore, one can easily show thatdet( D x f ( λ, − I ) = λ − c c ( b · a ) . Therefore, the Jacobian matrix of f ( λ, · ) at the origin has a simple eigenvalue of 1 if and onlyif λ equals λ ∗ = c c ( b · a ) . Now denote the right and left eigenvectors of D x f ( λ ∗ , x ∗ ) by v and w , respectively. One canshow directly that v = a defined in (6), and w is such that w = b · a , w d = b d , and w k = b k + S k w k +1 for k = 2 , . . . , d − . Then in order to establish the transcritical bifurcation, two nondegeneracy conditions haveto be verified. Since we have w t D λ f ( λ ∗ , x ∗ ) = 0, one first has to show that w t D xλ f ( λ ∗ , x ∗ ) v = ( b · a ) c c w t b is nonzero, which is clearly satisfied since all the terms of b and w are non-negative, andcontains terms of the form b k (which are strictly positive for each nonzero b k ).Second, we need to show that w t D xx f ( λ ∗ , x ∗ )[ v, v ] (cid:54) = 0. Since only the first componentof f , which we call f , is nonlinear, one merely needs to consider the second derivative of thiscomponent function. We get the following formula. D xx f ( λ ∗ , x ∗ )[ v, v ] = 2( β − α )Ω d (cid:88) k =2 p k a k . By looking at the corresponding parameter values, this value is also nonzero, and thereforethe second nondegeneracy condition holds. This completes the proof of the theorem.
In this paper, we have considered an age-structured population model for red coral popu-lations with a parameter of fitness. When the fitness increases sufficiently, a set of stableinvariant closed curves of oscillating orbits form, and these stable curves persist for largevalues of the fitness parameter. It is not surprising that for small fitness parameters, solu-tions limit to extinction, but we see that even for large fitness, populations become extremelyvulnerable, as they limit to oscillation spending long period of time near extinction.The coral population model has a curve of fixed points containing a Neimark-Sacker,saddle-node, and transcritical bifucation point. We develop new methods based on previouscomputer-assisted proof methods and use these methods to validate the branch of fixed points,and the three bifurcation points. 29 cknowledgments
We would like to thank Konstantin Mischaikow for pointing us to this coral population model.This research was partially supported by NSF grant DMS-1407087. In addition, E.S. andT.W. were partially supported by the Simons Foundation under Awards 636383 and 581334,respectively.
References [1] R. Aubourg. Red coral in the Mediterranean sea.
Wikimedia Commons , 2 November 2017.[2] S. Beslin, D. Baney, and V. de Angelis. Small denominators: No small problem.
Mathematics Magazine ,71(2):132–138, 1998.[3] L. Bramanti, M. Iannelli, and G. Santangelo. Mathematical modelling for conservation and managementof gorgonians corals: youngs and olds, could they coexist?
Ecological Modelling , 220(21):2851–2856, 2009.[4] L. Bramanti, G. Magagnini, L. D. Maio, and G. Santangelo. Recruitment, early survival and growth ofthe Mediterranean red coral Corallium rubrum (L 1758), a 4-year study.
Journal of Experimental MarineBiology and Ecology , 314(1):69–78, 2005.[5] M. J. Capinski, E. Fleurantin, and J. D. Mireles James. Computer assisted proofs of two-dimensionalattracting invariant tori for odes.
Discrete and Continuous Dynamical Systems, Series A , 2020. Toappear.[6] M. G. Crandall and P. H. Rabinowitz. Bifurcation from simple eigenvalues.
Journal of FunctionalAnalysis , 8:321–340, 1971.[7] S. Das, C. Dock, Y. Saiki, M. Salgado-Flores, E. Sander, J. Wu, and J. Yorke. Measuring quasiperiodicity.
Europhysics Letters , 114(4):40005, 2016.[8] P. G´ery. Corallium rubrum (Linnaeus, 1758) - Banyuls-sur-Mer, Sec de R´ed´eris: 08/84.
WikimediaCommons , 31 July 2011.[9] W. J. F. Govaerts.
Numerical methods for bifurcations of dynamical equilibria . Society for Industrial andApplied Mathematics (SIAM), Philadelphia, PA, 2000.[10] H. B. Keller.
Lectures on numerical methods in bifurcation problems , volume 79 of
Tata Institute ofFundamental Research Lectures on Mathematics and Physics . Published for the Tata Institute of Funda-mental Research, Bombay; by Springer-Verlag, Berlin, 1987. With notes by A. K. Nandakumaran andMythily Ramaswamy.[11] Y. A. Kuznetsov.
Elements of Applied Bifurcation Theory . Springer-Verlag, New York, second edition,1998.[12] J.-P. Lessard, E. Sander, and T. Wanner. Rigorous continuation of bifurcation points in the diblockcopolymer equation.
Journal of Computational Dynamics , 4(1–2):71–118, 2017.[13] S. M. Rump. INTLAB - INTerval LABoratory. In T. Csendes, editor,
Developments in Reliable Com-puting
Physica D , 411:132569, 2020.[15] E. Sander and T. Wanner. Validated saddle-node bifurcations and applications to lattice dynamicalsystems.
SIAM Journal on Applied Dynamical Systems , 15(3):1690–1733, 2016.[16] E. Sander and T. Wanner. Equilibrium validation in models for pattern formation based on Sobolevembeddings.
Discrete and Continuous Dynamical Systems, Series B , 2020. To appear.[17] G. Santangelo, L. Bramanti, and M. Iannelli. Population dynamics and conservation biology of theover-exploited Mediterranean red coral.
Journal of Theoretical Biology , 244(3):416–423, 2007.
18] J. B. van den Berg, J.-P. Lessard, and E. Queirolo. Rigorous verification of Hopf bifurcations via desin-gularization and continuation. arXiv:2006.13373 [math.DS] , 2020.[19] T. Wanner. Computer-assisted equilibrium validation for the diblock copolymer model.
Discrete andContinuous Dynamical Systems, Series A , 37(2):1075–1107, 2017.[20] T. Wanner. Computer-assisted bifurcation diagram validation and applications in materials science.
Proceedings of Symposia in Applied Mathematics , 74:123–174, 2018., 74:123–174, 2018.