Period-3 dominant phase synchronisation of Zelkova serrata: border-collision bifurcation observed in a plant population

The population synchrony of tree seed production has attracted widespread attention in agriculture, forestry and ecosystem management. Oaks usually show synchronisation of irregular or intermittent sequences of acorn production, which is termed ‘masting’. Tree crops such as citrus and pistachio show a clear two-year cycle (period-2) termed ‘alternate bearing’. We identified period-3 dominant phase synchronisation in a population of Zelkova serrata. As ‘period-3’ is known to provide evidence to imply chaos in nonlinear science, the observed period-3 phase synchronisation of Zelkova serrata is an attractive real-world phenomenon that warrants investigation in terms of nonlinear dynamics. Using the Hilbert transform, we proposed a procedure to determine the fractions of periods underlying the survey data and distinguished the on-year (high yield year) and the off-year (low yield year) of the masting. We quantified the effects of pollen coupling, common environmental noise and individual variability on the phase synchronisation and demonstrated how the period-3 synchronisation emerges through a border-collision bifurcation process. In this paper, we propose a model that can describe diverse behaviours of seed production observed in many different tree species by changing its parameters.

The highly synchronised fluctuation of annual seed production is common in perennial plant species. Among such species, acorn masting has particularly attracted interest in many fields of research [1][2][3][4] . In silviculture, prediction of masting behaviour in forest stands is necessary for approaches that promote natural regeneration 5,6 . As acorns are a substantial food source for wild animals 2 , it is important to understand masting for ecological management [7][8][9][10] . ' Alternate bearing' refers to tree crops that produce heavy crops one year (the 'on' year) and light crops the following year (the 'off ' year). Citrus (e.g., oranges, lemons and mandarins), pistachios and chestnuts are crops that show pronounced alternate bearing [11][12][13][14][15] . Acorn masting and alternate bearing also have been investigated in terms of the synchronisation of ensembles of trees [16][17][18][19] . In nonlinear physics, the synchronisation of ensembles of oscillators is known to be caused by mutual coupling or common identical noise 20,21 . Many types of coupling, such as indirect global and local coupling [16][17][18][19] and direct coupling 22 , have been investigated. The common noise-induced synchrony [23][24][25][26][27] is known as the Moran effect in population ecology. In a 15-year field survey, we observed period-3 dominant phase synchronisation in a population of 106 individuals of Zelkova serrata. To date, the majority of masting behaviour has been recognised as irregular and/or intermittent sequences and alternate bearing of tree crops is generally period-2. Therefore, the period-3 dominant synchronisation identified here is unique. In particular, period-3 is a special term in nonlinear dynamics, as it has been proven that certain dynamics of period-3 can generate any periods including chaos 28 . Thus, elucidation of the mechanism of period-3 phase synchronisation in Zelkova serrata will contribute to understanding the variety of periods synchrony reported in many perennial plant species.
The objective of this study was to clarify the mechanism underlying the period-3 dominant synchronisation in the population of Zelkova serrata surveyed. We developed a method to quantify various periodic compositions coexisting in an ensemble time series. We demonstrated a globally coupled map of the resource budget model (hereafter GCM-RBM) 16 to model the period-3 dominant phase synchronisation of Zelkova serrata as its masting behaviour.
The GCM-RBM has been commonly used to model the population synchrony of cross-pollinated plants [16][17][18] . At a certain magnitude of coupling strengths (β), the group consisting of N RBMs is strongly synchronised so that its dynamics can be represented by a one-dimensional map. We demonstrate how the period-3 emerges from a tangential bifurcation of the GCM-RBM whose map (characterised as a border-collision bifurcation) is piecewise smooth and piecewise monotonic 29,30 . We also estimated control parameters (R C and β) of the GCM-RBM for the survey data for Zelkova serrata. The proposed approach is expected to be a powerful method to understand the mechanisms of synchronisation in various perennial plant species.

Materials and Methods
Field survey. Zelkova serrata is a diclinous monoecious tree distributed in East Asia and flowers in April and May. The leaf colour changes to vivid red in autumn and the leaves on fruiting twigs change colour much earlier than those on non-fruiting twigs (see Supplementary Fig. S1). On the basis of visual inspection from the ground in two weeks of mid-November, the seed production level was classified into 10 classes. This method is popular in vegetation surveys 31,32 . Data from a population of 106 trees acquired over 15 years from 2003 to 2017 were analysed. This primary survey was conducted in the area bounded by longitude 139°28′45.57″E to 139°28′45.84″E and latitude 35°40′11.99″N to 35°40′31.90″N in Fuchu City, Tokyo, Japan. We also conducted a additional survey for 48 trees since 2006 in the area bounded by longitude 139°41′27.64″E to 139°41′40.58″E and latitude 35°41′25.04″N to 35°41′27.30″N in Shinjuku District located 20 km west of Fuchu City.
Fraction of the period-Q sequence. For identification of the 'on-year' and the 'off-year' for a single time series x(t), we defined the flag index ON(t) as a step function by employing phase-based and amplitude-based definitions.
As the phase angle θ(t) of a single time series x(t) is needed to define ON(t), it is given by where HT is the Hilbert transform of the true signal, X is the time average of X(t) [31][32][33][34][35][36][37][38] .
In the phase-based definition, we used the phase π 2 as a threshold to distinguish the 'on-year' and the 'off-year' . t ( ) 2 θ < π and θ ≥ π t ( ) 2 correspond to the on-year and the off-year, respectively.
In the amplitude-based definition, we used the time average of x(t) as a threshold to distinguish the 'on-year' and the 'off-year' . When x(t) is larger than , then year t is considered to be the 'on-year' .
A Matching between ON A (t) and ON P (t) was examined with 1600 combinations of the two control parameters (β and R C ) used in Fig. 7 to show that the two indices were identical.
The afore-mentioned two definitions of ON(t) can be expanded to the ensemble {x i (t); t = 1, 2, …, T, i = 1, 2, …, N} to obtain ON i (t) for each tree i as follows, www.nature.com/scientificreports www.nature.com/scientificreports/ To investigate the phase synchronisations of masting in a population of trees, it is important to observe the composition of periods in individual trees and the population. Here, we propose a practical means to determine the composition of periodic components. For instance, period-2 and period-3 sequences are defined as 'ON ⇒ OFF ⇒ ON' and 'ON ⇒ OFF ⇒ OFF ⇒ ON' , respectively. Thus, a period-Q sequence is defined as the sequence where one 'on-year' at year t is followed by Q − 1 for 'off-year' and 'on-year' arises at year t + Q. The fraction of period-Q in the i th tree's time series is determined by where ON i (t) is determined by Eqs (4,5). The fraction of period-Q representing a population is given by i The median and mode can be used to represent a population as well as the mean.
Measures of synchrony. Given that we focused on the phase synchronisation, we employed the notion of in-phase and out-of-phase analysis 22 . The fraction of in-phase in a population {x i (t); i = 1, 2, …, N, t = 1, 2…, T} is also used 22 to measure the phase synchronisation of a population of trees. If the arbitrary pair of x i (t) and x j (t) show in-phase behaviour between year t and year t + 1, then Let x i (t) be the yield of the i th tree in year t, and define φ (i, j, t) as the phase between the i th and j th tree: The fractions of in-phase movements of the i th tree relative to the remaining trees in the population is defined as where H IN is given as Heaviside step function, however, Model. The resource budget model (RBM) is described below and has been used previously to model masting and alternate bearing 16,39 . Let S i (t) be the amount of resource reserves at the beginning of year t for tree i. If the accumulated resource + S t P ( ) i S exceeds the threshold of the pool (L T ), then the excess amount S i (t) + P S − L T is used for flowering C f i .
The cost of pollinating flowers and bearing fruit is designated by C a i . The cost ratio After the reproductive stage, the accumulated resource becomes The RBM is a one-dimensional map modelled by Eqs (10)(11)(12), where −R C is the slope at the fixed point of the RBM. Isagi (1996) introduced the fruiting efficiency of a tree, Y(t) 16 , as a global coupling term.
where β is the strength of pollen coupling and N denotes the population size. Equation (12) is replaced by Eq. (14) to model the pollen coupling: a i C f i = As described above, the GCM-RBM is established with Eqs (10)(11)(12)(13)(14).
To model the phase synchronisation associated with a nontrivial disturbance, we incorporate individual noise (e I ) and common noise (e C ) into the GCM-RBM. Individual noise (e I ) assumes the heterogeneity of trees 24 . Common noise (e C ) can induce synchrony, which is known as the Moran effect 25 .
These noise types are imposed on P S in the following manner: Here, the random number σ(t) is drawn from the normal random number N (μ, σ 2 ) = (0, 1), and δ ii (t) is assigned to each tree individually. P 0 is the intrinsic annual surplus. The seed production level x i (t) obtained in the survey ( Fig. 1(a)) is considered proportional to in the GCM-RBM. In the numerical simulations, P 0 and L T were set as 10 and 100, respectively, (2019) 9:15568 | https://doi.org/10.1038/s41598-019-50815-8 www.nature.com/scientificreports www.nature.com/scientificreports/

Results and Discussion
Field experiments. Figure 1 shows the results obtained in the primary survey at Fuchu City. The seed production level of the individual trees, x i (t) and annual (ensemble) mean of the population at year t, are shown in Fig. 1(a). X (t) mostly shows period-3 cycles of an on-year (high-production year) followed by two consecutive off-years (low-production years). Figure 1 Fig. 1(e) in red, which shows the phase histogram for the 15-year term, clearly demonstrates the three evenly distributed peaks corresponding to the three fundamental phases (0, 2π/3 and −2π/3). Figure 1(f) displays a circle map in each panel in which blue dots are plotted at the points of (cosθ i (t), sinθ i (t)). On each panel, the arrow represents the vector (mean{cosθ i (t)}, mean{sinθ i (t)}), and the amplitude of the arrow is the order parameter that measures the strength of synchronisation. In the panel for 2003 in Fig. 1(f), the amplitude of the arrow is almost zero. This indicates that desynchronisation occurred in 2003. Therefore, this year should not be an 'on-year' even though Θ(2003) is almost zero. This is consistent with the median of {ON P,i (t), i = 1, 2, … N} indicating an 'off-year' in t = 2003. Hence, the phase θ(t) is useful to diagnose the states of the synchronisation of ensembles of trees. Thus, in Fig. 1, we detected a period-3 dominant phase synchronisation of Zelkova serrata. Period-2 synchronisation is common as alternate bearing in crop production, such as that of citrus [39][40][41] and nuts 42 . In addition, in acorn masting, irregular and/or intermittent sequences are also common 1,43,44 . Therefore, the dominance of the period-3 sequence observed in Zelkova serrata is unique and notable. The masting of Zelkova serrata is characterised by two key features: (a) two significant periods, i.e., period-3 and period-2 coexist, and (b) the fraction of period-3, FP (3), is significantly larger than that of period-2, FP (2). Figure 2 shows the results of the additional survey at Shinjuku District. The key features observed in the primary survey at Fuchu City shown in Fig. 1 are completely consistent with those of the additional survey at Shinjuku District. In particular, the agreements identified in the histograms of seed production (Figs 1(b) vs 2(b)) and the histograms of phase (Figs 1(b) vs 2(b)) are remarkable. The results indicate the presence of a long-range spatial synchronisation between the two populations at Fuchu City and Shinjuku District. The ON-OFF sequence of the two populations was perfectly matched from 2006 to 2017. The presence of such strong spatial synchronisation suggests that global pollen coupling occurred in the range of at least 20 km. The long-ranged spatial synchronisation of masting also has attracted widespread interest 45 . Figure 3 plots the density bifurcation diagrams of S and θ Cs for three β values. For β = 0, in Fig. 3(a), the first bifurcation occurs at R C = 1, and at the even integer values of R C , such as 2, 4, …, periodic solutions appear for S. For β = 6, in Fig. 3(b), the bifurcation diagrams of S display period-(Q + 1) windows in the interval where R C = Q belongs 17 . The period of the periodic window increases as R C increases. This bifurcation is www.nature.com/scientificreports www.nature.com/scientificreports/ known as the period-adding bifurcation, which is a typical property of border-collision bifurcation explained in Fig. 4 46,47 . The period-adding sequence is clearly demonstrated in Fig. 3(b,f). The density bifurcation diagrams of θ Cs are also displayed in Fig. 3(c,d) for β = 0 and 6, respectively. It is difficult to measure S for a real tree but it is possible to measure C S . Therefore, the phase (θ Cs ) is a powerful variable to investigate the periods of real-world masting data. The fractions of period-Q {FP(Q); Q = 2 …, 6} defined by Eq. (7) are expressed in Fig. 3(e,f) for β = 0 and 6, respectively. With pollen coupling of β = 6, FP(Q) is dominant (mostly 1.0) in the interval of the period-Q window, as shown in Fig. 3(f). This result indicates that the perfect period-(R C + 1) phase synchronisation arises in the ranges including every digit number of R C 17 . In the ranges between adjoining period-windows, we can determine the composition of various period-Q (Q = 2, 3 …) with FP(Q). The strength of phase synchronisation is estimated as F IN . For example, in Fig. 3(e), at up F IN is 0.5 or smaller because, there is no-coupling as β = 0. Contrarily, for β = 6, F IN is maintained at 1 for the entire range of R C shown in Fig. 3(f), which indicates perfect phase synchronisation.

Mechanism of period-3.
For β = 6, a clear period-3 window is present in the range 1.6172 ≤ R C ≤ 2.0 ( Fig. 3(b)). Therefore, we selected R C = 1.6171 to generate Fig. 4. The first iterated map is plotted in Fig. 4(a) with thick black dots for β = 6. The first iterated map for β = 0 is drawn in a thin blue line, which is a tent map as formulated by Eq. (10). The two first iterated maps of Fig. 4(a) are identical only in S(t) < L T . At the border, S(t) = L T , the left derivative is 1 for both β = 0 and 6. The right derivatives are −R C and 0 for β = 0 and 6, respectively. The third iterated maps for β = 0 and 6 are plotted in Fig. 4(b). The map for β = 6 exhibits a clear tangent bifurcation with the three tangency points of 84.4523, 94.4532 and 99.9438. These three tangency points are expressed with three circle marks. Figure 4(c,d) magnify the ranges where the two tangency points of 99.9438 and 84.4523 locate, respectively. The period-3 window of the GCM-RBM is explained by a tangent bifurcation 48,49 . Given that the map (β = 6) illustrated in Fig. 4(a) is a piecewise smooth and piecewise monotonic map, the bifurcation is the border-collision bifurcation causing the period-adding bifurcation 29,30 as observed in Fig. 3(b,f). This border-collision bifurcation has been studied in the dynamics of switching circuits 29,30 . The existence of the identically same dynamics in electronic circuits and perennial plants implies the universality of the boundary-collision bifurcation in nature.
Effects of noise on the GCM-RBM. The effect of the two types of noise (e C ) and (e I ) on the GCM-RBM was evaluated with four combinations of e C and e I at R C = 2 and β = 6 in Fig. 5(a-d). In the noise-free condition (Fig. 5(a)), perfect period-3 phase synchronisation is apparent. In the case of only common noise, e C = 0.2 ( Fig. 5(b)), perfect phase synchronisation is maintained, but four or more consecutive off-year sequences (i.e., larger periodic, irregular and/or intermittent) appear. It should be noted that imposition of common noise (e C ) does not generate any individual disturbances.   Figure 5(d) demonstrates significant disturbances; however, the appearance of four or more consecutive off-year sequences is inconsistent with the key features of the survey data for Zelkova serrata. It is obvious that the individual variability e I is essential, and the common noise e C is unnecessary, to explain the period-3 dominant synchronisation of Zelkova serrata. Conversely, irregular, intermittent sequences and/or long consecutive off-years have been reported in many other tree species, such as spruce, beech and oak 50 . These reported features correspond to the common noise-imposed case exhibited in Fig. 5(d). Thus, the imposition of both e I and e C could account for various masting behaviours 51,52 .
In Fig. 6, the density bifurcation diagrams with individual noise e I = 0.2 are plotted in the same arrangements as in Fig. 3. It is difficult to obtain clear structural information from S in Fig. 6(a,b) because of noise presence (e I = 0.2). However, even with noise, the density bifurcation diagrams of θ Cs show clear structural information (Fig. 6(c,d)) indicating the presence of three fundamental phases in the wide range of R C . By varying R C , three states are observed in Fig. 6(d,f). In state I, the majority of RBMs behave randomly such that F IN is much smaller than 1.0. In state II, F IN attains 1.0 and the period-2 cycle dominates. In state III, period-3 and larger periods emerge, and the fraction of period-2 FP(2) gradually declines. The densities around the three fundamental phases are significantly high in state III as shown in Fig. 6(d).
Parameter studies: estimation of β and R C for the field data. To estimate the range of R C and β for the survey data for Zelkova serrata, we conducted a parameter study and present the results in the β-R C diagrams of Fig. 7. First, we show the noise-free conditions in Fig. 7(a-c) corresponding to the fractions of period-2 FP(2), period-3 FP(3) and period-4 FP(4), respectively. The region of FP(2) = 1 is sickle-shaped, and those of FP(3) = 1 and FP(4) = 1 are wedge shaped.
For a noise-imposed condition (e I = 0.2), FP(2), FP (3) and FP (4) are plotted in Fig.7(d-f), respectively. F IN and FP(3)/FP (2) are shown in Fig. 7(g,h). The values FP(3) = 0.5354 and FP(3)/FP(2) = 2.879 were calculated for the primary survey data in Fig. 1. In Fig. 7(e), the range of FP(3) = 0.5354 is located above the sickle-shaped region. FP(3)/FP(2) = 2.879 also appears in the sickle-shaped region (Fig. 7(h)) and indicates the coexistence of period-3 and period-2 to the same extent as in Fig. 1. The strength of synchronisation is quantified by F IN . F IN is 1 in state II. The β-R C diagram of F IN is divided into three regions corresponding to the states I, II and III (Fig. 7(g)). State I is the desynchronised region. In state II, FP(2) and F IN are near or equal to 1.0, thus indicating almost perfect period-2 synchronisation. State III has phase synchronisation consisting of several different periods. The perfect period-(R C + 1) synchronisation is apparent in every digit number of R C 17 owing to the period-adding nature of border-collision bifurcation. With increasing noise level e I , the perfect phase synchronisation moves into imperfect phase synchronisations accompanied with partial desynchronisation. Some fractions of period-(R C + 1) are still located in the region around the digit numbers of R C and coexist with other periods. For example, the relatively higher intensity region of FP(2), FP (3) and FP(4) locates around R C = 1, 2 and 3, as shown in Fig. 7(d-f), respectively. As discussed above, we estimated that R C = 2, β = 6 are e I = 0.2 are appropriate values for Zelkova serrata. Figure 8 shows the results of a numerical experiment of the GCM-RBM with R C = 2, β = 6, e I = 0.2 and e C = 0. The first key feature of period-3 synchronisation, which is an on-year followed by two consecutive off-years, is clearly demonstrated as most of the seed production (C S ) behaves as a period-3 sequence (Fig. 8(a)). In Fig. 8(b), 70% of the counts were in the two lowest classes in the last panel. Fraction of periods, FP i (Q) and FP(Q), were www.nature.com/scientificreports www.nature.com/scientificreports/ demonstrated in Fig. 8(c) for individual trees and the population, respectively. The FP(Q) is determined as FP(2) = 0.260, FP(3) = 0.730, FP(4) = 0.020, FP(5) = 0.010 and FP(6) = 0.000. Period-3 is dominant followed by period-2 as FP(3)/FP(2) = 2.879, and the other periods are trivial. Three fundamental phases (0, 2π/3, and −2π/3) are repeated ( Fig. 8(d,e)). In the last subplots of Fig. 8(e), the phase θ iCs is evenly located around the three fundamental phases. Importantly, the results presented in Fig. 8 are consistent with those of Figs 1 and 2.
Concluding remarks. Synchronisation of seed production is ubiquitous in many tree species. In this paper, we clarified the mechanism underlying the period-3 dominant synchronisation in populations of Zelkova serrata surveyed in Tokyo. We developed a method to determine various periodic compositions coexisting in an ensemble time series. With this method, we found that the phase synchronisations of the populations of Zelkova serrata in Fuchu City and Shinjuku District are identical. The observed long-range spatial synchronisation implies global pollen coupling occurred in the range of at least 20 km.
We employed the GCM-RBM to describe the synchrony of the masting of Zelkova serrata. When the coupling of the GCM-RBM is sufficiently strong, the dynamics of the GCM-RBM become a piecewise smooth and piecewise monotonic map that is characterised by border-collision bifurcation and the period-adding bifurcation. The mechanism generating period-3 is explained as the tangential bifurcation of the GCM-RBM. The GCM-RBM realised the period-3 dominant phase synchronisation and the key features of the survey data.
The presence of a period-3 solution implies the coexistence of various periodic solutions and chaos in the system. In addition, the GCM-RBM generates periodic solutions larger than period-2 because of its period-adding bifurcation mechanism. Thus, the developed GCM-RBM shows the potential to describe diverse seed-production behaviours observed in many tree species by manipulating its control parameters (R C and β) and the levels of individual and common noise imposed (e I and e C ).