Universality of stable multi-cluster periodic solutions in a population model of the cell cycle with negative feedback

We study a population model where cells in one part of the cell cycle may affect the progress of cells in another part. If the influence, or feedback, from one part to another is negative, simulations of the model almost always result in multiple temporal clusters formed by groups of cells. We study regions in parameter space where periodic ‘k-cyclic’ solutions are stable. The regions of stability coincide with sub-triangles on which certain events occur in a fixed order. For boundary sub-triangles with order ‘ ’, we prove that the k-cyclic periodic solution is asymptotically stable if the index of the sub-triangle is relatively prime with respect to the number of clusters k and neutrally stable otherwise. For negative linear feedback, we prove that the interior of the parameter set is covered by stable sub-triangles, i.e. a stable k-cyclic solution always exists for some k. We observe numerically that the result also holds for many forms of nonlinear feedback, but may break down in extreme cases.


Background
In biology, there are phenomena that may emerge when cells or other dynamical units are coupled and influence each other. One example is synchronization, a state in which units are progressing with identical or nearly identical coordinates. In this manuscript, we consider a related phenomenon, phase synchronization or temporal clustering, in which multiple groups (or clusters) of cells are synchronized. They have the same phase as the other cells in their group, but the phase differs from one group to another. Throughout this manuscript, 'clustering' will mean temporal clustering (not spatial clustering).
There is a huge literature about synchronization, but temporal clustering has been studied much less often. It was observed in models of connected neurons [1,14,15,20,38] and models of coupled engineered biological oscillators [11,40]. Temporal clustering has been observed in some experiments, including coupled chemical oscillators based on the Belousov-Zhabotinski reaction [34,35] and engineered electrochemical arrays [16,17,37].
Temporal clustering also has been observed in Yeast Metabolic Oscillations (YMO) experiments that exhibit stable periodic oscillations of budding yeast (Saccharomyces cerevisiae) between aerobic metabolism and anaerobic metabolism. These have been observed in experiments and studied for many years [8,12,19,21,27,31].
The Cell Division Cycle (CDC) of budding yeast is the cyclic process of cell growth and division. A correlation between YMO and bud index, fraction of cells budded, was presented very early [19,21], but the link between YMO and the CDC was unclear due to the fact that the periods of YMO were always shorter than the CDC times in the same experiments. The relationship between YMO and the CDC was mostly ignored until it was noted in genetic expression data [18,31] and investigators began calling such oscillations 'cell cycle related'. Boczko et al. [3] proposed that temporal clustering is the mechanism linking the metabolism and the CDC in these experiments. Temporal clustering with two clusters was observed in YMO experiments reported in Stowers et al. [33] and Young et al. [39] and it was postulated that cells are segregated into CDC synchronized clusters due to feedback effects on the CDC progression. They proposed that cells in some phase of the cell cycle might produce signalling agents or metabolic products and the levels of these chemicals may affect the growth rate of cells in some other parts of the cycle.
Population level models of the cell cycle that include the effect of such feedback between different parts of the cell cycle were proposed in Boczko et al. [3] and Young et al. [39]. It was observed in [5] that clustering in these models is heavily related to a number theoretic relationship between the number of clusters k and other indices that characterize possible clustered solutions. The main goal of this manuscript is to prove a conjecture made in [5] that for negative feedback, temporal clustering is universal, meaning that it happens for any parameters in the model.

A cell cycle population model, notations and previous results
We will represent the progress of a cell by a variable (phase) z i ∈ [0, 1]; 0 is the beginning of the cell division cycle and 1 is the end of the cycle. We will identify 0 with 1 and consider z i as a normalized coordinate on the circle. We will let S = [0, s) and R = [r, 1) where 0 < s < r < 1; see Figure 1. The model of progress of the ith cell in the cell cycle conceptualized in [3] is as follows: where i = 1, . . . , n, I ≡ #{i : z i ∈ S} n (fraction of cells in the signalling region).
We will assume that f is a monotone function satisfying f (0) = 0 and f (I) > −1. It is perhaps nonlinear.
The key feature of the model is that cells in one part of the cycle (denoted by S = [0, s) for signalling) influence the cells in another part (R = [r, 1) for responsive). We postulate that . According to the model, z 1 (t) to z 7 (t) are progressing with rate 1 while z 8 (t) to z 11 (t) are progressing with rate 1 + f ( 3 11 ).
the feedback is not direct from cell to cell, but is exerted through the influence of chemicals in the medium that are produced or perhaps consumed by cells in different regions of the cycle. The numbers 0 < s < r < 1 are fixed for a given application. We will treat them as parameters of the dynamical system. Previous work [39] shows that the number of clusters in asymptotically stable solutions is strongly dependent on s and r.
We do not know how the coordinates of our model correspond to particular phases of the actual cell cycle and so it does not reflect specific biology. However, it is well established that different parts of the cell cycle produce, consume and react differently to different metabolites and other chemicals (see, e.g. [7,9,13,[24][25][26]36,41]). A checkpoint, a phase of the cell cycle where progress may be arrested until sufficient conditions are satisfied [13], is a natural candidate for a responsive phase (R in our model) with negative feedback. This particular mechanism was modelled and studied in [22].
We note that the right-hand side of the differential equation (1) is piecewise constant. Solutions fail to be differentiable at discrete points in time, and the definition of solution must be generalized to continuous, piece-wise smooth solutions. This was addressed in [39]. Figure 1 provides a visualization of the system in model (1) for the case of n = 11 where 11 different locations of cells represent the state of each cell in the cycle.
Note that cells are identical in (1) and so if two cells have the same coordinates at some time t , then they will coincide for all t > t . A cluster will mean a group of cells whose coordinates are equal and the integer k will be the number of clusters.
When the cells are grouped in k clusters it is sufficient to consider only the k coordinates, {X i } k i=1 of the clusters rather than the coordinates of all individual cells. It follows trivially that the progress of clusters is given by for i = 1, . . . , k, and I is again given by (2). If the k clusters all have the same number of cells, then the progress of evenly clustered solutions can be described by (3), with the fraction of clusters in the signalling region. Next, we define a special class of periodic, evenly clustered solutions of the system (3)-(4) that will be the main object of our study. Definition 1.1: Let n be the total number of cells in the cycle and k ≥ 2 be a divisor of n, with n/k cells in each cluster. Suppose that there exists a smallest positive number d such that X j (d) = X j+1 (0) for all j = 1, . . . , k − 1 and X k (d) = X 1 (0). Then {(X 1 (t), . . . , X k (t)) | t ∈ R + } is called a k-cyclic solution.
The existence of k-cyclic solutions was proved in [6,39]: For any f (I) > −1 and any pair (s, r), if k is a divisor of n, then a k-cyclic solution exists consisting of n/k cells in each cluster. For negative f and a given pair (s, r), the k-cyclic solution is unique, up to a translation in time.
The main result of this manuscript is that for negative linear feedback, given any pair of parameters (s, r), 0 < s < r < 1, there is a k ≥ 2 such that the unique k-cyclic solution corresponding to (s, r) is (locally) asymptotically stable. In order to prove this result, we will first establish new results about stability of k-cyclic solutions for (s, r) in certain subsets. This will require the introduction of quite a bit of terminology.
Let x i denote the initial condition of X i (t) and let 1 := {(x 1 , . . . , x k )|x 1 = 0}. For a solution with an initial condition in 1 , let t * be the shortest positive time required for X 1 (t) to return back to its starting position, 0; i.e. when the solution returns to 1 . The Poincaré map : 1 → 1 for the flow (3)-(4) is defined as follows: Because we have assumed that f (I) > −1, this Poincaré map is well defined.

Definition 1.2:
We define a map F by where t k is the shortest time required for X k (t) to reach 1.
It was noted that the map F k is conjugate to the Poincaré map [39]. Further, if (x * 2 , . . . , x * k ) is a fixed point of F, then (0, x * 2 , . . . , x * k ) ∈ 1 is the initial condition of a k-cyclic solution of (3)-(4) and vice versa. Moses [23] showed that for k even clusters, asymptotic stability in the k clustered phase space implies asymptotic stability in the full phase space of n cells, i.e. system (1)- (2). Thus, to study the stability of any k-cyclic solution, it is sufficient to study the stability of the fixed points of F.
Note that 0 = x 1 ≤ x 2 ≤ · · · ≤ x k ≤ 1 can be considered as a (closed) simplex. The map F is defined on the interior of the simplex, i.e. on 0 < x 2 < x 3 < · · · < x k < 1. It was shown in [39] that the map F can be extended continuously to the boundary of the simplex and that this continuously extended F permutes components of the boundary of the simplex.
Define σ to be the number of clusters in S = [0, s) at the initial time and define ρ to be the number of clusters in [0, r) (outside of R) at the initial time. That is, According to Equations (3)-(4), clusters may change their rate only when one of them reaches s, r, or 1 which we call events s, r and 1, respectively.
be a solution whose initial condition satisfies (7). Then, X σ (t) reaching s is called the event s, X ρ (t) reaching r is called the event r, X k (t) reaching 1 is called the event1. A sequence sr1 . . . means that the clusters progress with an order of events s, r and 1, respectively. If two or three events happen at the same time, we say that the solution has simultaneous events.
The sequence of events followed by a cyclic solution is periodic and either sr1sr1 . . . or rs1rs1 . . . or has simultaneous events [5]. Thus, if we wish to study the stability of cyclic solutions, we need only to consider two periodic sequences of events, sr1 or rs1.
We will call the set of points {(s, r) : 0 < s < r < 1} the parameter triangle .

Definition 1.4:
A subset τ ⊂ is called isosequential for k if the k-cyclic solution corresponding to each parameter pair in the interior of τ has the same σ and ρ and the same order of events: That is, s (σ , ρ, k) and r (σ , ρ, k) are regions on which the k-cyclic solutions have the order of events sr1 or rs1, respectively. Figure 2 shows indices of isosequential regions for the case of k = 3. Isosequential regions are in fact sub-triangles [5]. Figure 3 illustrates that the relative positions of the parameters and the clusters of the cyclic solution determine the order of events.
Isosequential regions are important because it was shown in [5] that all the k-cyclic solutions for (s, r) in one of these sub-triangles have the same stability type. We will call a sub-triangle stable if all the k-cyclic solutions corresponding to points in the region are asymptotically stable. An isosequential region will be called unstable if all the k-cyclic solutions corresponding to points in the region are unstable. An isosequential region is called neutral if all the k-cyclic solutions corresponding to points in the region are neutrally stable (stable, but not asymptotically stable). In the context of k-cyclic solutions, stable means that initial conditions near the orbit of the periodic solution will remain near the orbit for all forward time. Asymptotic stability means that nearby solutions not only remain in a neighbourhood of the periodic orbit but also converge to the periodic orbit in forward time.
A point (s, r) ∈ is called a simultaneous point if the k-cyclic solution has events s, r and 1 occurring simultaneously. Note that if (s, r) is a simultaneous point, then the initial  Top: X 1 (t) will reach s before X 2 (t) reaches r, so the pair (s, r) ∈ s (1,2,3). Bottom: X 2 (t) will reach r before X 1 (t) reaches s, so (s, r) ∈ r (1,2,3).
condition (assumed to be in 1 ) of the k-cyclic solution has one cluster located at s and another one at r. For any f and integer k ≥ 2 isosequential regions are sub-triangles with simultaneous points at their corners [5]. In the interior of each sub-triangle, all k-cyclic solutions have the same order of events and same stability.   [5]. The left panel corresponds to a positive feedback: the red sub-triangles were shown to be unstable in [2,23], the green regions are unstable for large feedback [23]. We will show that they are neutral for small feedback. The right panel corresponds to a negative feedback: the blue regions were shown to be asymptotically stable and the yellow areas neutrally stable in [23], the red regions are unstable [2].
By the boundary sub-triangles, we mean the vertical, horizontal and oblique sub-triangles. Figure 4 shows sub-triangles according to this definition for the case of k = 5.
Next we review previous results about the stability of cyclic solutions. The vertical and horizontal boundary sub-triangles with order of events sr1, i.e. s (σ , k, k) and s (1, ρ, k), were shown to be neutral for any f (I) > −1 in [5].
For positive feedback, all interior sub-triangles are unstable and for negative feedback all interior sub-triangles with the order sr1 are unstable [2,23].
In conclusion, stability under negative and positive feedbacks of the vertical and horizontal sub-triangles was fully investigated in earlier works. However, stability of the oblique sub-triangles has not been previously analysed. The oblique sub-triangles are indexed by All results that had been established previously are shown in Figure 5 for the case of k = 8 where stability of the white subtriangles has not been investigated before. In this manuscript, we prove that the stability of the oblique triangles follows a similar pattern as the vertical and horizontal sub-triangles.

Main results
As discussed above, the stability of the vertical and horizontal sub-triangles was completely investigated, but the stability of the oblique sub-triangles was open. The ideas presented in Section 3 complete the study of stability of boundary sub-triangles. The stable sub-triangles are those with order of events rs1 and whose index is relatively prime with respect to the number of clusters k. Results in Section 3 can be combined into the following which was conjectured in [5]. This result is illustrated in Figure 6 for k = 8. The grey (sr1) and white (rs1) areas in Figure 6 are neutral. The blue shaded areas (rs1) in Figure 6 are stable.
Our method not only clarifies the stability of oblique sub-triangles but also provides an alternative proof for the vertical and horizontal sub-triangles. We establish these results by analysing all eigenvalues of the Jacobian matrix of the map F corresponding to points in each isosequential region.
In Section 5, we prove our main result, which is the following.
Theorem B: For any negative linear feedback f (I) = −cI, c < 1, and any pair (s, r) in the interior of the triangle 0 < s < r < 1, there is a positive integer k, k ≥ 2, such that the k-cyclic solution corresponding to (s, r) is stable. This result is important because it means that negative linear feedback systems of this type will always possess an asymptotically stable clustered periodic solution, irrespective of the values of parameters.
This was conjectured more generally in [5]. Numerically, this result appears to hold also for non-linear negative feedback, except possibly for some extreme cases. (See the last section for discussion.) In Section 4. we prove a version of Theorem thmB in the case of the zero feedback limit where the main ideas are more clear.
These results appeared in expanded form in the dissertation of the first author [28].

Some preliminaries
We introduce the following notations that will be used for the rest of this manuscript.

Order of events rs1
Suppose that σ < ρ and let {(X 1 (t), . . . , X k (t)) | t ∈ R + } be a solution of (3)-(4) such that its initial condition (where X i (0) = x i ) satisfies (7) and This condition implies that the solution has order of events rs1.
The determinant of DF rs1 was calculated to be which is less than one in absolute value since −1 < ω < 0. For 1 ≤ σ < ρ ≤ k, eigenvalues of DF rs1 for sub-triangles r (σ , ρ, k) satisfy the equation and 1 is not an eigenvalue [2].

Order of events sr1
Consider a k-clustered solution with order of events sr1. Then, its initial condition satisfies At the time t k , [2] obtained the solution This defines the map F sr1 . The Jacobian matrix of the map DF for the order of events sr1 was shown in [2] to be the following: We note that this implies that k cyclic solutions with order sr1 cannot be asymptotically stable. For 1 ≤ σ ≤ ρ ≤ k, eigenvalues of DF sr1 for sub-triangles s (σ , ρ, k) satisfy the equation and 1 is not an eigenvalue [2].  By the symmetry of zeros on the unit disc, each element ofA i has a complex conjugate in Recall that for a complex number, z = x + iy, an argument of z, arg(z), is a number φ ∈ R such that x = r cos(φ) and y = r sin(φ), where r = |z|. Figure 7 shows the general picture of arg(z − z 1 ), arg(z − z 2 ), and arg( z−z 2 z−z 1 ) for given z, z 1 , z 2 ∈ C. The argument of z − z 1 is the angle between the line from z to z 1 and the horizontal line passing z 1 .
The following lemma might be a known fact, but since we cannot find it in the literature, we state it here.
(depending on the location of the arc A). If λ / ∈ A, (depending on the location of the arc A).
We will use the following theorem found in [4] repeatedly.

Theorem 2.2 ([4]):
All roots of the polynomial p(s) = a n s n + a n−1 s n−1 + · · · + a 1 s + a 0 , a n = 0 and a i ∈ R for i ∈ {1, . . . , n}, are in the interior of the unit ball if and only if • | a 0 a n | < 1, The following lemma will help us when we use Theorem 2.2.

Lemma 2.3:
Let p 1 (λ) and p 2 (λ) be real monic polynomials of even degree n with their zeros contained in S 1 . Let a 1 , . . . , a n be zeros of p 1 (λ) and b 1 , . . . , b n be zeros of p 2 (λ), where a i = b j (1 ≤ i, j ≤ n) and Let A i be the smaller open arc of S 1 bounded by a i and b i for i = 1, . . . , n. Let θ be a real number such that θ ∈ (−1, 0) ∪ (0, 1). If the arcs A i are pairwise disjoint and each element of A i has a complex conjugate in A n−i+1 (i = 1, . . . , n), then there is exactly one zero of in each A i for any positive real number k 0 .

Proof:
The key idea of this proof is found in [10].
First, let θ ∈ (0, 1) and k 0 ∈ R + . Assume that the arcs A i are disjoint and each element of A i has a complex conjugate in A n−i+1 (i = 1, . . . , n). Define a function For fixed i ∈ {1, . . . , n}, we claim that A(λ) is a continuous real-valued function on each arc A i . Let λ ∈ A i . Consider Let θ i be the angle of the arc A i . By Lemma 2.1, Thus arg(q i q n−i+1 ) = π . This implies q i q n−i+1 is a negative real number. For j ∈ {1, . . . , n} and j = i, Lemma 2.1 implies Thus arg(q j q n−j+1 ) = 0. This implies q j q n−j+1 is a positive real number. Hence, is a negative real number, so for any λ ∈ A i . Then, A(λ) is a continuous real-valued function of λ on the arc A i . Since A(λ) has the values 0 and 1 at the endpoints of the arc A i , A(λ) must take all values between 0 and 1 on the arc A i . By the Intermediate Value Theorem, there exists z ∈ A i such that This implies that so z is a zero of (12). This means that there exists a zero of (12) in each A i for i = 1, . . . , n.
Since there are n different zeros of (12), there is exactly one zero of (12) in each A i (i = 1, . . . , n).

Classification of sub-triangles for negative feedback
For a negative feedback system, we have −1 < ω < 0. For the reader's convenience, we provide a table of symbols in Appendix.
We now consider the stability of oblique sub-triangles.

Theorem 3.1:
For negative feedback and k ≥ 2, If gcd(k, σ ) = 1, then the sub-triangles Proof: By substituting ρ = σ + 1 in Equation (9), eigenvalues of DF rs1 for sub-triangle and 1 is not an eigenvalue. Let θ = −ω. In a negative feedback system, we get 0 < θ < 1. Let Our claim is that all roots of f θ (λ) satisfy |λ| < 1. We will use Theorem 2.2 to obtain the claim. First, multiply both sides of (14) by (λ − 1) 2 , Note that f θ (λ) is equal to a polynomial of degree k−1, except at λ = 1. Since we consider λ = 1 in (15), we have f θ (λ) = f θ (λ), so we can refer to both as f θ (λ). The constant term of f θ (λ) is 1 and the leading coefficient of Since | 1−θ 1 | < 1, we achieve the first condition of Theorem 2.2. Let Then Next, we compute D 1 (λ, f θ ) and D 2 (λ, f θ ). Since we get Hence, Similarly, To check the second condition of the theorem, we consider two cases of the integer k: it is an odd number and it is an even number. Case I: k is an odd number.
In this case, k − σ can be either an odd or even number. Note that regardless of whether k − σ is an odd number or an even number, −1 is a zero of f 1 (λ). If k − σ is an odd number, then −1 is a zero of f 1 (λ) since σ is must be an even number, see (16). If k − σ is an even number, −1 is a zero of f 1 (λ), see (16). Since gcd(σ , k − σ ) = 1 and k is an odd number, we can write We use the following lemma which we will prove later.

Lemma 3.2:
For every 1 ≤ i ≤ k − 2, the arc C i contains at most one zero of f 1 (λ).
Let A i be the smaller open arc of S 1 bounded by Note that b k−1 2 = −1. By our construction, A i are disjoint and each element of A i has a complex conjugate in A k−i for i = 1, . . . , k − 1 (see Figure 8). Lemma 2.3 implies has exactly one zero in each A i (i = 1, . . . , k − 1). Hence, all zeros of D 1 (λ, f θ ) lie on S 1 . By the construction of A i (see Figure 8), the zeros of D 1 (λ, f θ ) alternate with the zeros of D 2 (λ, f θ ) for all θ ∈ (0, 1). Therefore, by Theorem 2.2, all zeros of f θ (λ) are in the interior of the unit disc, |λ| < 1, for all θ ∈ (0, 1). Case II: k is an even number. In this case, σ and k − σ are odd numbers and Let Note that p 1 (λ) and p 2 (λ) have no common zeros and they both have the same degree k−2, which is an even number. We write Note that b i = −1, 1 and c j = −1, 1. Let C i be the smaller arc of S 1 as follows: Note that p 1 (λ) has k−2 2 roots whose arguments are between 0 and π . By Lemma 3.2, the arcs C i each contain one root of has exactly one root in each A i (i = 1, . . . , k − 1). Hence, all zeros of D 1 (λ, f θ ) lie on S 1 . By the construction of A i in this case, the zeros of D 1 (λ, f θ ) alternate with the zeros of D 2 (λ, f θ ) for all θ ∈ (0, 1). Therefore, by Theorem 2.2, all zeros of f θ (λ) satisfy |λ| < 1, for all θ ∈ (0, 1).
The following theorem was proved in [23] with a different method.
The next corollary is proved by similar ideas as in the proofs of Theorem 3.1. In this corollary, there are some eigenvalues of DF rs1 lying on the unit disc, |λ| = 1.

Corollary 3.4:
For negative feedback with the order of events rs1 and k ≥ 2, the following statements hold.
Suppose that gcd(k − σ , σ ) = n = 1. We consider (15) in the proof of Theorem 3.1, We claim that all zeros of g θ (λ) lie in the unit disc. When the claim holds, it follows that no zero of f θ (λ) is outside the units disc, i.e. |λ| ≤ 1. Therefore, r (σ , σ + 1, k) are neutrally stable. Since gcd(k, σ , k − σ ) = n, the function g θ (λ) is a polynomial function of degree k−n and g θ (λ) satisfies the first condition of Theorem 2.2. Let Then, Next, we compute D 1 (λ, g θ ) and D 2 (λ, g θ ) defined in Theorem 2.2. Since Hence, Similarly, To check the second condition of Theorem 2.2, we consider three cases: k is an odd number, k is an even number and n is an odd number, and, k and n are both even. Case I: k is an odd number. This case implies that n is an odd number. Then k−n is an even number. Since σ and k − σ cannot be both odd numbers in this case, −1 is a zero of g 1 (λ). By our construction, g 1 (λ) and g 2 (λ) have no common zero. We can write Let A i be the smaller open arc of S 1 bounded by Note that b k−n 2 = −1. By our construction, A i are disjoint and each element of A i has a complex conjugate in A k−n−i+1 (i = 1, . . . , k − n). Lemma 2.3 implies that has exactly one zero in each A i (i = 1, . . . , k − 1). Hence, all zeros of D 1 (λ, g θ ) lie on S 1 and the zeros of D 1 (λ, g θ ) alternate with the zeros of D 2 (λ, g θ ) for all θ ∈ (0, 1). Therefore, all zeros of g θ (λ) are in the interior of the unit disc, |λ| < 1, for all θ ∈ (0, 1). Case II: k is an even number and n is an odd number. In this case, σ and k − σ are odd numbers, k−n is an even number, and we have , Note that p 1 (λ) and p 2 (λ) have no common zeros and they both have the same degree k−n−1, which is an even number. We write has exactly one zero in each A i . Hence, all zeros of D 1 (λ, g θ ) lie on S 1 and the zeros of D 1 (λ, g θ ) alternate with the zeros of D 2 (λ, g θ ) for all θ ∈ (0, 1). By Theorem 2.2, all zeros of g θ (λ) are in the interior of the unit disc for all θ ∈ (0, 1). Case III: k and n are even numbers.
In this case, σ , k − σ , and k−n are even numbers, and we have Note that p 1 (λ) and p 2 (λ) have all zeros on S 1 and they have no common zeros. Both p 1 (λ) and p 2 (λ) have the same degree k−n, which is an even number. We write has exactly one zero in each A i (i = 1, . . . , k − n). Hence, all zeros of D 1 (λ, f θ ) lie on S 1 and the zeros of D 1 (λ, g θ ) alternate with the zeros of D 2 (λ, g θ ) for all θ ∈ (0, 1). Therefore, all zeros of g θ (λ) are in the interior of the unit disc for all θ ∈ (0, 1).
We prove the following lemma to investigate stability of the boundary sub-triangles corresponding to order of events sr1. Lemma 3.5: Let k and σ be positive integer satisfying 1 ≤ σ ≤ k and let n = gcd(k, k − σ + 1) and m = gcd(k, σ ), where k + 1−n−m is an even number. Let If p 1 (λ) and p 2 (λ) have no common zeros and −1, 1 are not their zeros, then the function has all zeros on the unit disc for all θ ∈ (0, 1) ∪ (−1, 0). Proof: If we substitute σ = 1 or ρ = k in (11), the left-hand side of (11) becomes zero. Thus the eigenvalues of DF sr1 for the regions s (1, σ , k) and s (σ , k, k) satisfy
• k is even and n is odd. Since n = gcd(k, k − σ + 1) and n is odd, we have that k − σ + 1 is odd. Then, σ is even. Since m = gcd(k, σ ), it follows that m is even. In this case, k, σ , m are even numbers and n, k − σ + 1 are odd numbers. • k and n are even. Since gcd(n, m) = 1 and n is even, it follows that m is odd. Then, σ is odd due to m = gcd(k, σ ) and k is even. We then get that k − σ + 1 is odd. In this case, k and n are even numbers and m, σ , k − σ + 1 are odd numbers. Let Case I: k, n, m, σ and k − σ + 1 are odd numbers.
In this case, we obtain that k + 1−n−m is an even number. Moreover, p 1 (λ) and p 2 (λ) have no common zeros and −1, 1 are not their zeros. By Lemma 3.5, all zeros of P θ (λ) lie on S 1 .
Case II: k, σ , m are even numbers and n, k − σ + 1 are odd numbers. By Lemma 3.5, all zeros of P θ (λ) lie on S 1 . Case III: k and n are even numbers and m, σ , k − σ + 1 are odd numbers. By Lemma 3.5, all zeros of P θ (λ) lie on S 1 .
From this proof, we also get the following corollary: Corollary 3.7: Let k and σ be positive integers satisfying 1 ≤ σ ≤ k and let θ ∈ (0, 1) ∪ (−1, 0). The equation has all solutions on the unit disc.
The results in this section identify the location of all asymptotically stable and neutrally stable regions for the boundary sub-triangles. See Figure 10 for an illustration.
We illustrate stability regions of the boundary sub-triangles under negative feedback in Figures 11 and 12. For a given feedback function f (I) and a given k, a program computes vertices of each boundary sub-triangle and uses the results in the previous section to check stability of each sub-triangle. The vertices of each sub-triangle can be computed by a formula presented in the next section.    Figure 11 shows the stability of boundary sub-triangles of for the case of k = 10 with different negative feedback functions.
In Figure 12, we show stability regions for the cases of k = 23 (prime) and k = 24 (composite) under the different negative feedback functions.
We see in these figures that for different feedback functions, the pattern of stability persists, while the sub-triangles become skewed.

The universality of cyclic solutions
In the previous section, we identified which boundary sub-triangles are (asymptotically) stable and which are neutral. In this section, we study how the sub-triangles fit together. For a given linear negative feedback, we prove that the interior of the is covered by the asymptotically stable sub-triangles.
Previously, we saw that an isosequential region has a triangular shape inside and its vertices are simultaneous points [5].
We first find the formula of each simultaneous point. A simultaneous point c σ ,ρ,k is a point (s, r) in such that the k-cyclic solution corresponding to the point satisfies the following: • it has events s, r, and 1 occurring simultaneously, • the solution at the initial time has σ clusters in the S region and it has ρ clusters outside the R region. Figure 13 illustrates the location of clusters at corner points of r (σ , ρ, k). Let {(X 1 (t), . . . , X k (t)) | t ∈ R + } be the k-cyclic solution that corresponds to c σ ,ρ,k . We will continue to use x i = X i (0). Then, the k-cyclic solution at the initial time satisfies that r = x ρ+1 and s = x σ +1 . The clusters outside the R region are equally distributed by a distance d (as in Definition 1.1 of k-cyclic solution) and the clusters inside the R region are equally distributed by a distance called d , see Figure 14. It follows that ρd + (k − p)d = 1. Note that distances outside of R correspond directly to time. Let t ρ be the time required for X ρ (t) to reach r = x ρ+1 and t k be the time for X k (t) to reach 1. We will seek conditions for t ρ = t k . Recall that β σ ,k = f (σ/k). Thus If we require t ρ = t k and since ρd + (k − p)d = 1, we have We now consider the case of no feedback, i.e. f (I) = 0 for all I. We include this section in order to illustrate some of the main ideas of the proof of Theorem B in a context that is more clear.   With zero feedback, events and isosequential regions are still defined. The formula of each simultaneous point, (34), for f ≡ 0 becomes simply Definition 4.1: Let k, σ and ρ be positive integers such that 1 ≤ σ < k and 2 ≤ ρ ≤ k. A sub-triangle r (1, ρ, k) is a relatively prime sub-triangle if gcd(ρ − 1, k) = 1. If gcd(σ , k) = 1, sub-triangles r (σ , k, k) and r (σ , σ + 1, k) are also called relatively prime sub-triangles.
For nonzero negative feedback, relatively prime sub-triangles are the ones that are stable as shown in Theorems 3.1 and 3.3.
The next lemma is the main idea of the proof. It states that if a boundary triangle is not relatively stable, then it is included entirely inside a larger boundary triangle that is relatively prime. This idea is illustrated in Figure 16. The overlay of relatively prime sub-triangles for k = 2 (yellow), 3 (red), 4 (green) and 6 (black) over k = 12 (blue). The key feature is that these relatively prime sub-triangles perfectly cover the white rs1 boundary sub-triangles on the left.

Lemma 4.3: Consider a zero feedback system. Let (s, r) be a point in the parameter triangle . Then there exist a positive integer k ≥ 3 and k such that the point (s, r) is in the triangle k .
Note here that the following statements hold for 2 ≤ i ≤ k − 1: • c 1,i,k are on the line r = 1 k , The next lemma is the second important idea of the proof. We use induction on k to show that for each k, the 'gap' k+1 \ k is covered by relatively prime subtriangles. Specifically, the gap is covered by boundary triangles of the form r (σ , ρ, k) and r (σ , ρ, k − 1). While not all of these triangles are relatively prime, they are all included perfectly in a larger relatively prime sub-triangle by Lemma 4.2.

Lemma 4.4: Consider a zero feedback system. For a positive integer k ≥ 3, k is covered by relatively prime sub-triangles.
Proof: We prove this lemma by induction on k. For k = 3, 3 = {c 1,2,3 } ⊆ (1, 2, 2) which is a relatively prime sub-triangle. Suppose that k is covered by relatively prime sub-triangles. It suffices to find a relatively prime covering of the gap between k+1 and k , that is the region k+1 \ k . The gap can be considered as composed of the vertical, horizontal and diagonal gaps. First, we consider the vertical gap between k+1 and k . Let j be a positive integer such that 1 ≤ j ≤ k − 2. Let A be the point on the lines s = 1 k . We will show below that ABC is covered by the triangle r (1, j + 1, k − 1). See Figure 18 for a visualization of ABC. Note that, By the definition of points A, B and C, we get Moreover, the following inequalities hold: This inequality implies that the r-coordinate of A is greater than or equal j/(k − 1). Since it follows that the slope of line passing j/(k − 1) and B is less than or equal to 1. Hence, (39), (40), and (41) imply that ABC is covered by the sub-triangle r (1, j + 1, k − 1), see Figure 19. Since r (1, j + 1, k − 1) is covered by relatively prime sub-triangle (see (36) in Lemma 4.2), the entire vertical gap between k+1 and k is covered by relatively prime sub-triangles. A similar idea can be used for the horizontal gap. Next, we consider the oblique gap, the area inside bounded above by r = s + 1/k and bounded below by r = s + 1/(k + 1), see Figure 20. Let j be a positive integer such that 1 ≤ j ≤ k − 1. Let a be the point on the lines s = j k and r = s + 1 k and b be the point on  k . We claim that abc is covered by triangle r (j, j + 1, k − 1). See Figure 20 as a visualization of abc. Note that We get the formulas for the point a, b, c as follows: Inequality (42) implies that the s-coordinate of c is less than or equal to j/(k − 1). By inequality (40), we have that the r-coordinate of a is greater than or equal to j/(k − 1). Hence, abc in Figure 20 is covered by triangle r (j, j + 1, k − 1).
Combining Lemmas 4.3 and 4.4, we get the following main theorem.

Covering under negative linear feedback
In this section, treat negative linear feedback function, f (x) = −cx for 0 < c ≤ 1.
In the previous section, we have proved under the zero feedback that the parameter triangle is covered by relatively prime sub-triangles. The main idea was that boundary sub-triangles that are not relatively prime are included perfectly in larger relatively prime sub-triangles. This fact was used to show that the gap between k+1 and k is covered by a union of relatively prime sub-triangles.
The calculation here is more complicated since the right angle triangles in Section 4 become a scalene triangles under negative linear feedback and Equations (36) and (37) in Lemma 4.2 no longer hold. Only the oblique boundary sub-triangles still have this property (Theorem 5.5). For vertical and horizontal sub-triangles, there are thin slices that are not covered by the same larger triangles as in the zero feedback case. The first technical difficulty is showing that these slices are in fact covered by other relatively prime sub-triangles (Lemma 5.9). The second difficulty is showing that these skewed sub-triangles fit together nicely to cover the gap, which is bounded by skewed lines (Lemma 5.10). Using Equation (34), the formula for a simultaneous point c σ ,ρ,k is as follows: From this formula, we obtain the following four propositions.
This proposition provides the fact that all simultaneous points c σ ,i,k lie on the same line for fixed σ and i = σ , . . . , k, see Figure 22 (the black line).

Proposition 5.2: Consider a negative linear feedback system. Let ρ and k be positive integers
This proposition provides the fact that all simultaneous points c i,ρ,k lie on the same line for fixed ρ and i = 0, . . . , ρ, see Figure 22 (the blue line).  The next result clarifies that any oblique sub-triangle corresponding to the order of events rs1 is covered by stable sub-triangles.
Note that c 0,ρ,k = c 0, ρ n , k n . By (43), Then m 1 ≥ m 2 . Hence, r (1, ρ + 1, k) ∩ k+1 is partially covered by stable sub-triangle r (1, p n + 1, k n ). Specifically, r (1, ρ + 1, k) ∩ k+1 can be considered as having two distinct parts; the part that intersects r (1, p n + 1, k n ) (the blue shaded area in Figure 27) and the part that intersects r (1, ρ + 1, k)\ r (1, p n + 1, k n ) (the red shaded region in Figure 27). The blue shaded area is covered by the stable triangle r (1, p n + 1, k n ). We claim that the red shaded area is covered by another stable sub-triangle. There are two cases depending on ρ to show that the red shaded area in Figure 27 is covered by a stable sub-triangle.
Case I: 2 ≤ ρ ≤ 2 3 (k − 1). Consider sub-triangle r (1, ρ m + 1, k−1 m ) where gcd(k − 1, ρ) = m. We claim that the red shaded area is a subset of r (1, ρ m + 1, k−1 m ).  Let m 6 , m 7 , m 8 be the slopes of l 6 , l 7 , l 8 , respectively, see Figure 28. Let A be the crossing point between l 2 and l 5 and let B be the crossing point between l 5 and l 8 , see Figure 28. To get the claim, we need to show that m 8 < m 7 < m 6 and the r-coordinate of A is greater than the r-coordinate of B. By (43), the slopes of m 6 , m 7 and m 8 are as follows: It follows that m 8 ≤ m 7 ≤ m 6 . The point A is the solution of the systems Then, the r-coordinate of A is The point B is the solution of the systems Then, the r-coordinate of B is We claim that r A − r B > 0. By simplification, it is equivalent to show that implies the inequality (45). Note that n = gcd(ρ, k) ≤ k 2 , so n − 2 3 ≤ k 2 − 2 3 . Then, implies the previous inequality. The left-hand side of the inequality (47) can be considered as a quadratic polynomial with variable c and fixed k. Then, the following properties hold: • coefficient of c 2 is nonnegative, • slope of the polynomial is negative at c = 0 and 1, • the values of the function are positive at c = 0 and 1.
Note that the left boundary of each r ( σ −1 l , k−2 l , k−2 l ) and r ( σ n , k n , k n ) has slope less than 1 by Proposition 5.4. To prove the claim that the blue region in Figure 30 is a subset of , it suffices to compare the s-coordinates of points A 1 , A 3 and A 4 , i.e. show that Figure 30 is a visualization of the claim. By the definitions of the points A 1 and A 3 , we get It follows that s A 1 < s A 3 for 2 ≤ σ ≤ k/5. Since we get Case II: k/5 ≤ σ ≤ k − 1. In this case, we claim that the blue region in Figure 30 is covered by r ( σ −1 q , k−1 q , k−1 q ). Let B 4 be the crossing point between r − k k+1 = ck (k+1) 2 s and the line passing through c σ −1,k−1,k−1 and c σ −1 . We claim that the s-coordinate of point A 3 is less than the s-coordinate of point B 4 . Figure 31 is a visualization of this claim. By Proposition 5.1 and the definition of B 4 , we get Since A < 0, it follows that This inequality implies that s A 3 < s B 3 . Therefore, the blue regions covered by stable subtriangles.

Lemma 5.9: Let f be a negative linear feedback function, f
The following will confirm that k (defined in Definition 5.6) is covered by stable subtriangles.   It suffices for the inductive step to show that the gap between k and k+1 , see Figure 33, is covered by stable sub-triangles. We separate the gap between k and k+1 into three parts as follows: • The vertical gap is the area bounded by four lines: the line passing through c 1,k+1,k+1 and c 1,1,k+1 , the line passing through c 1,k,k and c 1,1,k , the line passing through c 0,k,k+1 and c k,k,k+1 , and the line passing through c 0,1,k+1 and c k,k+1,k+1 . • The horizontal gap is the area bounded by four lines: the line passing through c 0,k,k+1 and c k,k,k+1 , the line passing through c 0,k−1,k and c k−1,k−1,k , the line passing through c 1,k+1,k+1 and c 1,1,k+1 , and the line passing through c k,k+1,k+1 and c 0,1,k+1 . • The oblique gap is the area bounded by four lines: the line passing through c 0,1,k and c k−1,k,k , the line passing through c 0,1,k+1 and c k,k+1,k+1 , the line passing through c 0,k,k+1 and c k,k,k+1 , and the line passing through c 1,k+1,k+1 and c 1,1,k+1 .
First, we consider the vertical gap and let 1 ≤ j ≤ k − 2.Let l = gcd(j, k − 1). Let l 1 , l 2 , l 3 and l 4 be lines as follows: • l 1 is the line passing through c 1,k+1,k+1 and c 1,1,k+1 , • l 2 is the line passing through c 0,j+1,k and c 1,j+1,k , • l 3 is the line passing through c 0,j,k and c 1,j+1,k , • l 4 is the line passing through c 0,j,k−1 and c 1, j l , k−1 l . By Lemma 5.9, the regions r (1, j + 2, k) ∩ k+1 and r (1, j + 1, k) ∩ k+1 (the blue shaded areas in Figure 34) are covered by stable sub-triangles. To show that the entire vertical gap is covered by stable sub-triangles, we further need to show that the region bounded by l 1 , l 2 and l 3 (the red area in Figure 34) is covered by r (1, j l + 1, k−1 l ). Note that r (1, j l + 1, k−1 l ) is a stable sub-triangle by Theorem 3.3. Let B 1 be the crossing point between l 1 and l 2 , A 1 be the crossing point between l 1 and l 3 and A 2 be the crossing point between l 1 and l 4 . To prove that the red area in Figure 34 is covered by r (1, j l + 1, k−1 l ), we show the following: Claim 1: The slope of the line passing through c 0,j,k−1 and B 1 is less than or equal 1.

Claim 2:
The r-coordinate of A 1 is greater than the r-coordinate of A 2 , i.e. r A 1 > r A 2 .
We first prove Claim 1. By (43), it follows that B 1 is the solution of Then, the r-coordinate of B 1 is crossing point between r − k k+1 = ck (k+1) 2 s and l 6 , and b σ be the crossing point between r − k k+1 = ck (k+1) 2 s and l 7 .Since the slope of the line passing through c σ where s a σ , s A σ , s B σ , s b σ denote the s-coordinates of a σ , A σ , B σ and b σ , respectively. Figure 35 is a visualization of this claim. Point a σ is the solution of Then, Then, Since B σ is the solution of Thus Lemma 5.12 implies s a σ < s A σ and s B σ < s b σ . Hence, the red shaded triangle in Figure 35 is covered by r ( σ q , k−1 q , k−1 q ) which is a stable sub-triangle. Therefore, the entire horizontal gap between k and k+1 is covered by stable sub-triangles.
Next, we consider the oblique gap between k and k+1 . There are four steps to prove how the oblique gap is covered by stable sub-triangles.
Step 1: We will show that the intersection of the vertical and oblique gaps and the intersection of the horizontal and oblique gaps are covered by r (1, 2, k) and r (k − 1, k, k). Let A 1 be the lowest corner of the former intersection and A 2 be the highest corner of the latter intersection. The two red shaded areas in Figure 36 are a visualization of the intersections. We note that the slope of the line segment from c 0,1,k to c 1,2,k is greater than 1. To complete this step, we need to show that A 1 ∈ r (1, 2, k) and A 2 ∈ r (k − 1, k, k). Point A 1 is the solution of r = s + 1 k + 1 and Then, Point A 2 is the solution of

Then,
, Let B 1 be the crossing point of r = s + 1 k+1 and the line passing through c 0,1,k and c 1,1,k and let B 2 be the crossing point of r = s + 1 k+1 and the line passing through c k−1,k,k and c k−1,k−1,k . To show A 1 ∈ r (1, 2, k) and A 2 ∈ r (k − 1, k, k), we need to show that ) .
Since Figure 37. A visualization of the claim to show that A 1 ∈ r (1, 2, k) and A 2 ∈ r (k − 1, k, k).
we get s A 1 > s B 1 . Since For Steps 2, 3 and 4 below, we define the following: • G σ denotes the region bounded by r = s + 1 Step 2: Let σ = 1, 2, k − 3, k − 2. The green regions shown in Figure 36 are a visualization of G 1 , G 2 , G k−3 , G k−2 . In this step, we claim the following: and We first prove ( . We claim that the points D σ and E σ are in r (σ , σ + 1, k − 2). Figure 38 is a visualization of this claim. We will prove this by comparing the s-coordinates of the points d σ , D σ and Then, The s-coordinate of d σ is as follows: By the formulas for s d σ and s D σ , we get s d σ < s D σ for σ = 1, 2. Point E σ is the solution of Then, We get the s-coordinate of e σ as follows: It is elementary to check that s E σ < s e σ for σ = 1, 2. Then, (59) holds. Next, we will prove (60). Let σ = k − 3, k − 2. Let d σ be the solution of r = s + 1 s and let e σ be the solution of r = s + 1 k+1 and . The s-coordinates of the points d σ , D σ , E σ , e σ are as follows: , , , Hence, we get (60).
Step 4: Let 3 ≤ σ ≤ k − 4 and k ≥ 11. In this step, we the proof will depend on the value of c. There are two cases to be considered; 0 ≤ c ≤ 0.5 and 0.5 ≤ c ≤ 1.
Case I: 0 ≤ c ≤ 0.5. We claim that the points D σ , E σ are in r (σ , σ + 1, k − 1). Let d σ be the solution of r = s + 1 k+1 and r − σ k−1 = cσ (k−1−σ ) (k−1) 2 s and let e σ be the solution of r = s + 1 k+1 and r − 1 = − (k−1−cσ )(k−1) Figure 40 is a visualization of   Point E σ is the solution of Then, The s-coordinate of e σ is as follows: By Lemma 5.14 below, we obtain s d σ < s D σ and s E σ < s e σ . Then, By Steps 1-4, the entire oblique gap between k and k+1 is covered by (67) Lemma 5.12: Let k and σ be positive integers such that 2 ≤ σ ≤ k − 3. Then, the following inequalities hold: for 0 < c ≤ 1.
The expression A can be rewritten as a quadratic polynomial in variable c for fixed k and σ : The polynomial has the following properties: • the coefficient of c 2 is negative, • its values at c = 0 and c = 1 are both negative, • the slopes at c = 0 and c = 1 are both positive.
Thus A < 0. It implies that Then, (68) holds. Next, we claim that Let For fixed k and σ , B can be rewritten as a quadratic polynomial in variable c: The polynomial has the following properties: • the coefficient of c 2 is positive, • its values at c = 0 and c = 1 are both negative, • the slopes at c = 0 and c = 1 are both positive.
Lemma 5.13: Let k and σ be positive integers such that k ≥ 11 and 3 ≤ σ ≤ k − 4. For 0 ≤ c ≤ 0.5, the following inequalities hold: Proof: Let g, h, i be functions on [3, k − 4] such that Note that By using the conditions of k and c, it is obvious that g (3) Next, let
Let a point (s, r) ∈ T. Then (s, r) ∈ k for some k. By Lemma 5.10, it follows that the point (s, r) is covered by a stable sub-triangle. We have now proved Theorem B.

Conclusion and discussion
In Theorem thmA, we have fully characterized the stability k-cyclic solutions for (s.r) in all boundary sub-triangles for any negative feedback. For the order of events sr1 they are always neutrally stable. For the order rs1, the stability is completely determined by the index i = σ or k − ρ. Solutions are asymptotically stable if i and k are relatively prime and neutral otherwise.
We have shown in Theorem thmB that for linear negative feedback the stable boundary triangles completely cover the parameter triangle .
Thus the model predicts that there will always exist an asymptotically stable clustered solution regardless of parameter values. Numerics show that this is also true for many forms of nonlinear feedback. See Figures 44 and 45. In [39], it was shown that for any negative feedback the synchronized solution is unstable. In [2], it was shown that the 'uniform solution', i.e. a steady-state periodic solution with cells spread out maximally around the circle is also unstable in most cases. Thus it appears that for negative feedback systems of this form and any parameter values, a k cluster solution will be the only stable periodic solution. This is what we always observe in simulations starting from random initial conditions -a k cluster periodic solution emerges with roughly the same number of cells in each cluster (see [5] and [32]). The universality of this result is important because for these models, as with many biological systems, the parameter values are hard to estimate from first principles.
In the following figures, we show visualizations for overlays of stable sub-triangles. The figures were generated using a MATLAB program by the following process: • We first consider linear negative feedback functions and then some examples of nonlinear negative feedback functions. • For each k where k = 2, . . . , m, the program constructs stable sub-triangles (blue shaded sub-triangles) by using the vertex formula (34) and Theorems 3.1 and 3.3. Note that these two results provide all stable sub-triangles among the boundary sub-triangles. • Then, the program displays in one plot all the constructed stable boundary subtriangles.    We provide two types of overlays; k ranges from 2 to 10 and k ranges from 2 to 100. Note here that the lighter blue regions in all figures throughout this section are the area that some parts of stable sub-triangles do not overlap other stable sub-triangles. These are the exception rather than the rule, indicating that for most parameter values there exist multiple stable cyclic solutions (with different values of k). In Figures 44 and 45, we show overlays of stable regions for two nonlinear feedback functions f (I) = −I/(2 − I) and f (I) = −3I 2 + 2I 3 . For both of these examples and many others that we examined, the parameter triangle is fully covered by stable sub-triangles. Thus we believe that the conclusion of Theorem thmB holds far beyond the linear case.
We note from Figures 16 and 21 and those in this section that most of the area of is already covered with stable sub-triangles with small k. While larger k are theoretically possible, they are expected to be of little consequence for three reasons: (1) the parameter sets where they exist grow small as k increases (they are on the boundary of ), (2) these sets are often overlapped by stability regions with fewer clusters and (3) the amount of feedback exerted by a single cluster would decrease like 1/k as k increases and at some point would be insufficient to overcome the noise in the system. The number of cells, n, in a bioreactor with a litre of fluid is on the order of 10 10 , so for practical purposes k n. We observe that there are some negative nonlinear feedback functions such that some points in are not covered by any stable boundary sub-triangle. By Theorems 3.1 and 3.3 and the vertex formula of sub-triangles (34), if there exist a positive integer k ≥ 2 and a subset G of k ( k is defined in Definition 5.6) such that G is not covered by stable boundary sub-triangles for l ≤ k, then G is not covered by any stable boundary sub-triangle. For f (I) = − √ I, Figure 46 suggests that there are regions in that are not covered by any stable boundary sub-triangles. This feedback function has the features that it is not Lipschitz and it is not bounded below by the functions covered in Theorem thmB. We also observe in Figure 47 that for a particular negative feedback function f (I) = −0.1 arctan(50(I − 0.22)) − 0.14801 that is Lipschitz and bounded below by the function −I, the union of stable boundary sub-triangles for all k ≥ 2 does not cover the interior of . Whether Theorem thmB can be extended more generally is an open question. The numerical results shown in Figures 46 and 47 make it clear that for some non-linear negative feedback functions the stable boundary sub-triangles alone do not fully cover the parameter triangle. As discussed in [2], there are a few exceptional interior subtriangles that are stable. It is possible that these stable interior sub-triangles cover any holes. In [29,30], the stability of exceptional triangles was proved for k = 9 and k = 14. Further, they showed that for k prime, some interior sub-triangles can become stable for larger feedback. The indices of these sub-triangles satisfy certain number theoretic relations. Since the holes in the figures appear only for larger feedback, it is possible that these exceptional sub-triangles cover some or all of the holes.
Numerical simulations, i.e. running the model for some (s, r) in the holes that appear show that even in the holes in the figures, the solutions always converge to some clustered solution.
In light of these considerations, settling the conjecture fully appears to be complicated.

Disclosure statement
No potential conflict of interest was reported by the author(s).