Dynamic Transitions and Stability for the Acetabularia Whorl Formation

Dynamical transitions of the Acetabularia whorl formation caused by outside calcium concentration is carefully analyzed using a chemical reaction diffusion model on a thin annulus. Restricting ourselves with Turing instabilities, we found all three types of transition, continuous, catastrophic and random can occur under different parameter regimes. Detailed linear analysis and numerical investigations are also provided. The main tool used in the transition analysis is Ma \&Wang's dynamical transition theory including the center manifold reduction.

1. Introduction. Acetabularia acetabulum is an important species for studying the morphogenesis during development for a single cell organism. In this article we will focus on the vegetative whorl generation along the stalk. Biological details of the Acetabularia development can be found in various review articles, for example [2] [3].
Several different models have been put forward to explain the whorl generation phenomenon. We note here since Acetabularia is a single cell organism and evidence has been shown that without nucleus the stem is still capable of growing a cap [4], the whorl generation is largely a local dynamics involving chemical reactions and cytoskeleton hence hope is given for a simple model with differential equations not involving genetics or intracellular signaling. Chemical as well as mechanical models [12][6] [7] have been proposed.
In this article we will apply the model used by Murray [13] , which is a purely chemical model and a reaction diffusion equation involving two reactants, and the 2 YIQIU MAO, DONGMING YAN AND CHUNHSIEN LU pattern formation is driven by the Turing instability of the system. Several observations have been shown to support the kinematic nature of our model: The spacing of the hair is affected by temperature [8], the spatial pattern of whorl hair is related to outside calcium concentration [5], and the whorl pattern is initiated simultaneously.
We will focus ourselves here with the effect of calcium concentration on the initiation of the whorl pattern. Experiments have shown [4] that there is a threshold of calcium concentration of about 2mM for the formation of whorls. Interestingly for external calcium above about 60mM, whorls will stop forming once again. Both facts can be shown in our model and details of the transition from stable stem to whorl formation when Calcium concentration crosses a critical value are explicitly calculated, using dynamical transition theory developed by Ma and Wang [10] [11].
More specifically, we obtain a specific interval of calcium concentration for the whorl pattern to occur. Then we proceed to classify the onset of pattern formation into three types, corresponding to the three types in Ma & Wang's theory. We have shown that all three kinds of transition type can occur though in the case of thin stem wall only continuous type and catastrophic type can happen. We also show in the numerical section that the number of whorl hairs generated during transition can be quite irregular with underlying parameters.
Besides the purpose of investigating the Acetabularia whorl formation, this paper can be viewed as an example for exploring dynamic behavior for reaction diffusion equations in an annular region.
The article is organized as following: the modeling background is introduced in Section 2, along with the deviation equation and functional setup and a discussion of global existence; in Section 3 we analyze the linear problem, from properties of Bessel functions we derive the crossing eigenvectors must be two dimensional when stem walls are thin, exact classifications of parameter space are provided, along with the principle of exchange of stabilities; in Section 4 we focus on the dynamical transitions of the system and prove three types of transitions can occur in different cases; in Section 5 we use numerical tools to investigate the transition types and properties of critical eigenvectors; in Section 6 we derive some physical conclusions from our research.
2. Model equations. In this paper, we apply a model proposed by Murray [13], which is an adaptation of a simple two species mechanism, the Schnackenberg (1979) system. The equations state as follows: Where all k's and D A , D B are all positive parameters, A and B are functions of r, θ (or x) and t with the annulus domain defined by A, B define the density of two substances inside the annular growth region at the tip of mature Acetabularia. Here we assumed a reaction 2A + B k3 − → 3A took place. In our context, B is calcium, and A could be a molecule that is increased in the calcium presence but are constantly transformed to other substances at a rate of k 2 . And they are generated by constant rate k 1 , k 4 respectively for A and B. We use the non-dimensionalization as follows: Then we have the reaction diffusion system after omitting the stars: where the Laplacian The δ is assumed to be very close to 1 as the stalks of Acetabularia have very thin walls.
The inward flux constant λ is directly proportional to outside calcium concentration because the normal intracellular concentration of calcium is extremely low compared to outside. Hence we will use λ as a principle parameter and investigate the transition when the increase of λ causes instability of the steady state.
It is easy to see that is a steady state of (2). Consider deviation of (2) from (3). To do so, we take Dropping the primes, the system (2) is then transformed into Next, we convert the equation (5) into an abstract functional setting which is standard in the framework of dynamic transitions. First, we let X and X 1 be function spaces defined as follows: Then, we define the operators L λ = −A + B λ : X 1 → X and G : X 1/2 → X (X 1/2 to be explained later) by where U = (u, v).
We thus obtain the equivalent operator equation of the evolution equation (5) as follows dU dt = L λ U + G(U ), It is easy to see that L λ is a completely continuous field, meaning it's the sum of a linear homeomorphism A : X 1 → X (with graph norm on X 1 ) and a compact operator B λ : X 1 → X, which insures relatively nice spectral property for L λ , see [10] Chap 3.
Since H 1 (Ω) 2 = X 1/2 is continuously imbedded in L 6 , thus it's also easy to check from (7) that the cubic G is locally Lipschitz continuous from X 1/2 to X which guarantees the local existence. Moreover G(U ) = o(||U || X 1/2 ) and is the sum of several continuous multilinear operators from X 1 to X. Hence the requirement on P527 of [11] is satisfied.
In [16], Y. You proved a global attractor (hence global existence) exists for this equation for any initial value U 0 ∈ X with Dirichlet boundary conditions.
For the second part, let's prove the following claims: 1. If we denote the smallest positive zeros of J n , J n , Y n , Y n as j n , j n , y n , y n , then we have n < j n < y n . 2. J n and Y n have alternating zeros which go to infinity when x > n. Proof of claim 1: [15] page 485-487 has shown that n < j n < y n < j n and J 2 n + Y 2 n is decreasing on the interval 0 < x < j n . Hence holds, using Y n (n) > 0 and Y n < 0 for n < x < j n , we derive that Y n (j n ) > 0 hence n < j n < y n . Proof of Claim 2: We know J n and Y n have alternating zeros that go to infinity, in between each consecutive zeros of J n or Y n there must be a zero for J n or Y n hence all the zeros must go to infinity. Furthermore it's straightforward to verify that Y n and J n satififies: Since Y n and J n has no common zeros (a property from linear independence between (Y n , Y n ) and (J n , J n ), (Y n , Y n ) and (J n , J n ) are independent as well. Hence by the Sturm-Liouville theory and claim 1, they have alternating zeros and all zeros are greater than n. Proof of lemma: We investigate the following quantity: This shows, along with claim 1, that the function Y n (x) first increase and peak at x = n and then decreases to −∞ at y n . And from y n on, by claim 2, this function decrease from +∞ to −∞ at each interval between consecutive zeros of Y n . Now the equation (13) can be rewritten as it then follows naturally that when δ is close to 1, λ n,1 lies in (n 2 /δ 2 , n 2 ) for n ≥ 1 while λ 0,1 = 0, hence λ n,1 → n 2 when δ → 1. While (δ − 1) λ n,2 is always greater than the minimum of distance between any two different zeros of Y n and J n that are smaller than δ λ n,2 , this forces λ n,2 to go to ∞ as δ → 1, and the case for any λ n,j , j > 1 follows true by the first part of the lemma. Remark 1. The above hypothesis resembles that of Bourget's hypothesis regarding zeros of Bessel function of the first kind, which was proved to be true [14][1]. Here the difference is our hypothesis focuses on the zeros of cross products involving the derivatives of the first and the second Bessel functions.
Given the above lemma and hypothesis, we can see that the multiplicity of λ k in (12) is two when n ≥ 1, and is one if n = 0. The hypothesis is not needed in the limiting case when δ is close to one.
Then the eigenvalue problem (9) can be solved using the ansatz Substituting (15) into the equation (9) , we obtain that From this we can see that β satisfies: Then Now we proceed to classify the parameters R, d, a, λ according to what unstable modes this linear equation generates while omitting the critical cases. We obtain a Proposition as follows: In what follows, we assume R, d, a, λ are all positive numbers. Here e k,l is stable (unstable) means the corresponding eigenvalue β k,l has negative (positive) real parts.
1. If (λ − a) > (λ + a) 3 , then the constant eigenvector e 0,l is unstable. 2 , then e k,2 is unstable, where λ k lies between: and all the other eigenvectors (including e 0,l ) are stable. If no λ k lies in this interval then every eigenvector is stable. + λ), then all the eigenvectors are stable.
Proof. Case 1: The equation (17) can be written as: It's easy to see if λ k = 0 and (λ − a) > (λ + a) 3 the equation has either two positive roots or two complex roots with positive real parts hence the conclusion. Case 2: The discriminant for the quadratic function viewed as a function of λ k is hence the second inequality appeared in case 2 guarantees the positiveness of the discriminant, since it also guarantees the coefficient for λ k in the quadratic function to be negative, basic quadratic analysis shows the for λ k between the two (positive) roots appeared in case 2, the constant term in (21) is negative hence guarantees a positive real root for (21). Case 3: The second inequality in case 3 shows the discriminant above is negative hence the constant term in (21) is always positive hence no positive real part is possible for any of the modes.
The parameter space in the case 2 of the above theorem can be shown between the two surfaces (outside the small cylinder and inside the bigger cone) in figure 1. Case 1 indicates the region inside the cylinder, while case 3 indicates the region outside both the cylinder and the ucone.
By treating λ as a changing parameter, which represents the concentration of calcium ions, we can expect that when λ increases from region 3 to region 2 in the above theorem, we should have a dynamic transition from a uniform steady state to a spatially variant steady state. We are particularly interested in this case because the other case where region 3 transitioned into region 1 is locally equivalent to an ODE version of (9) without the diffusion terms, since all spatial variant modes are stable during the transition. Actually we can prove:  Proof. This is a direct consequence of the Theorem 6.1.4 of [9]. With our invariant manifold being the constant function space.

YIQIU MAO, DONGMING YAN AND CHUNHSIEN LU
For simplicity, we assume a > 1 3 √ 3 . It's easy to verify under this condition, (λ − a) < (λ + a) 3 always hold for any λ > 0. Another simplification we will make is to restrict our transition purely inside region 2 instead of on the boundary of region 2. The reason behind this is the length of window for λ k is close to 0 near boundary and it is not generic that there is a λ k lies precisely inside this small window. The third simplification we shall make is there is only one (not counting multiplicity) λ k entering the length window when λ increase to a critical value. It's easy to see this is the generic case. The double entering case corresponds to boundary lines in Figure 2 later in numerical section hence is negligible.
Remark 2. Notice in the above theorem the multiplicity of β kc,2 depends on the underlying λ kc . If the corresponding n c (defined to be the n in the critical eigenvector R(r) cos nθ or R(r) sin nθ) is 0 then the multiplicity is 1, if n c > 0 then the multiplicity is 2. When δ → 1, according to Lemma 3.1 and boundedness of I λ when λ satisfies the restrictions in Theorem 3.4, only multiplicity 2 transition will occur.

Dynamical transition theorems. Suppose in this section the eigenvectors can be represented as
for n > 1, j ∈ Z + , l = 1 or 2, and for n = 0, and all the other quantities are denoted in the same fashion. (for example β njl to denote the eigenvalues). Here we assume all the eigenvectors are normalized, i.e. Ω |e sub | 2 = 1 for any subscripts, and e nj1c = e nj2c , e nj1c = e nj2c and e 0j1 = e 0j2 for complex eigenvectors, otherwise all eigenvectors are real. Moreover we can standardize the choice of u njl and v njl so that they have the following properties given the restrictions in the PES: The only non-trivial part of the above properties is the second one, which can be derived from (16), namely, Hence if we denote p l = v njl /u njl , then it's easy to derive p 1 · p 2 = 2λ (a + λ) 3 > 0, then the properties follows. See Theorem 2.1.3 of [11] for a classification of dynamical transition into three types.
Here we first define the following quantity q(λ): Recall n c is defined to be the n in the critical eigenvector R(r) cos nθ or R(r) sin nθ.
Notice that B nj1 = B nj2 for non-real eigenvectors hence q(λ) is a real number.
where Φ is the local center manifold function from space {y c e c + y s e s } to its complement in X. This equation will describe the same dynamic as the center manifold locally around 0. According to Theorem A.1.1 of [11], this center manifold can be approximated as follows: where P 2 is the standard projection onto the complement space, G 2 is the quadratic part of G, and o(2) = o(|y c | 2 + |y s | 2 ) + O(|β ncjc2 |(|y c | 2 + |y s | 2 )).
To get a third order approximation of (25), we only need second order terms in Φ, which is equal to Φ 2 (y c , y s ) = −L −1 λ P 2 G 2 (y c , y s ) by the above equation. It will be clear later that the only useful modes in Φ 2 are of follows: 1. Modes e 0jl for j ∈ Z + and l = 1, 2 and 2. Modes e 2ncjlc and e 2ncjls for j ∈ Z + and l = 1, 2. Using similiar calculation as above we obtain: and Now we go back to the equation (25), calculation shows Ω G(y c e c + y s e s + Φ(y c , y s )) · e c = Ω γ 2 [R ncjc (y c cos n c θ + y s sin n c θ)] 2 u 2 ncjc2 + γ 12 [R ncjc (y c cos n c θ + y s sin n c θ)] 2 u ncjc2 v ncjc2 + γ 3 [R ncjc (y c cos n c θ + y s sin n c θ)] 3 u 2 ncjc2 v ncjc2 + 2γ 2 [R ncjc (y c cos n c θ + y s sin n c θ)]u ncjc2 Φ u (y c , y s ) + γ 12 [R ncjc (y c cos n c θ + y s sin n c θ)]u ncjc2 Φ v (y c , y s ) where Φ u and Φ v are u and v components of Φ 2 , and o(3) is same notation as o(2). It's clear from this equation that only those modes with 2n c and 0 in Φ 2 wavenumbers are nonzero after integrating. After plug in with (26) where as is in the statement of the theorem. Hence the equation (25) can be rewritten as (by symmetry between y c and y s ) at β ncjc2 = 0, o(3) reduce to o(|y c | 3 + |y s | 3 ) and energy estimate gives Since u ncjc2 − v ncjc2 is always positive by previous discussion, the steady state 0 is locally asymptotically stable at critical λ if q(λ c ) is negative, which by Theorem 2.2.11 of [11] indicates a type I (continuous) type transition for the original system. If q(λ c ) > 0, then by Theorem 2.4.16 of [11], the system undergoes a type II (catastrophic) transition. 1. Equation (8) has a type III (random) transition from (0, λ c ). More precisely, there exists a neighborhood U ⊂ X of 0 such that U is separated into two disjoint open sets U λ 1 and U λ 2 by the stable manifold Γ λ of u = 0 such that the local transition structure is as shown below and the following hold: 2 to a unique singular point v λ on λ > λ c that is an attractor such that for every φ ∈ U λ 2 , lim t→∞ ||u(t, φ) − v λ || = 0.

Equation
Proof. Here we assume the transitioning mode to be then the projection of the center manifold is dy dt = β 0jc2 y + Ω G(ye + Φ(y)) · e, and .
then the results follows directly from Theorem 2.3.2 of [11].
5. Numerical Investigations. As is indicated in Remark 2, the transition when crossing eigenvalues are of multiplicity 2 is the most probable case when δ is close to 1, we are more interested in Theorem 4.1. However the coefficients p(λ c ) have to be computed numerically.
In this section we first show four plots of n c as a function of d and a, with R = 4, 8 and δ = 1.2 and 2. Then we calculate p(λ c ) for several selected points, and then show the graph for several critical eigenvectors. All coding are done through Mathematica. We are not varying R too much in our n c graph because the effect of R on critical eigenvectors are quite clear from the formula of I λ in Theorem 3.4. From the above graph we can see that critical wavenumber n c can change quite erratically when the parameters are near the boundary of the always stable regime. An interesting interpretation is the number of whorl hairs, which can be represented roughly as the number of positive patches of critical eigenfunctions, also behave irregularly. For example when R = 16, δ = 2 d = 80, if we increase a from 0.2 to 0.55, the number of whorl hairs during the initial whorl formation goes through 1, 2, 3,4,8,5,9,6,10,7,1,11,2,3,8,12,4,5,13. This shows the importance of distribution of Bessel function zeros since critical eigenvectors are affected by them.
The following is a chart for p(λ c ) for various parameters when the critical wavenumber is not zero. Majority of the parameters we have tested show jump transitions. Figure 3 shows a plot corresponding to a typical critical eigenvector with n c = 6, j c = 3.

6.
Conclusions. We discuss the effect of δ, R, λ, a, d on the system.   Figure 3. Critical Eigenvector cos n c θR nc,Jc (r) with n c = 6, j c = 3 and δ = 1.2, dark or bright regions indicate potential whorl hair growth, in this graph we have around 18 growth regions From the Remark 2, we can see that realistically, when the stem wall of Acetabularia is thin (δ is small), the crossing modes are two-dimensional and hence Theorem 4.1 applies and there is either a jump or continuous type of transition depends on q(λ c ). If the stem wall is thick, then a random type of transition can occur and the concentration for substance u or v can display a concentric pattern.
R affect the general reaction and diffusion rate for both substances, temperature can be a good candidate. From (22) and Figure 2, we can see that greater R (or reaction rate) corresponds to more whorl hair formation because higher wavenumber eigenvectors are involved. This correspond well to the experiment observation from Harrison [8].
The effect of λ as we are principlely concerned with, is causing a Turing type of transition of the system from a uniform steady state to a spatially variant steady state like those of Figure 3. However, since both continuous and catastrophic type transitions are found from our numerical investigation, the critical eigenvector only faithfully represent the transition state when the type is continuous. For the catastrophic type, we only know the solution will go away from the origin to a global attractor, the precise bifurcated solution is unknown. It is also worth mentioning that λ need to be located in a finite interval for the steady state to lose it's stability, as can be seen clearly from Figure 1.
Parameter a shows the generating rate for substance u. Our result shows that when this rate is high, the whole system will be stable and no transition occurs.
The amplitude of d is inversely proportional to the diffusion rate for u, by Theorem 3.2, d need to be > 1 for Turing instability to occur hence u should be a larger molecule than Calcium. Interestingly, our Figure 2 shows that the larger molecule u is, the fewer whorl hair will be generated during transition.