Estimation of Blooming Start with the Adaptation of the Uniﬁed Model for Three Apricot Cultivars ( Prunus armeniaca L.) Based on Long-Term Observations in Hungary (1994–2020)

: The aim of our research was to adapt Chuine’s uniﬁed model to estimate the beginning of blooming of three apricot cultivars (‘Cegl é di b í borkajszi’, ‘Gönci magyar kajszi’, and ‘R ó zsakajszi C.1406’) in Hungary in the time period 1994–2020. The uniﬁed model is based on the collection of chilling and forcing units. The complexity of the model lies in the high number of parameters necessary to run it. Following the work of other researchers, we reduced the number of relevant model parameters (MP) to six. In order to estimate the six MPs, we used a simulated annealing optimization method (known for being effective in avoiding getting stuck in local minima). From the results, we determined the local optimum of six MPs, and the global optimum parameter vector for three apricot cultivars. With these global optimum parameter vectors, the beginning of blooming could be estimated with a root-mean-square error (RMSE) of less than 2.5 days, using the knowledge of the daily mean temperature in the time period 1994–2020.


Introduction
The bud burst and the blooming are the most extensively studied phenophases of temperate trees (for more details see [1]). Such studies found that the dominant environmental factor affecting the growth speed of trees is air temperature [2]. In the last few decades, several temperature-based, mechanistic phenology models were produced for simulating the vegetative (e.g., bud burst) or reproductive (e.g., blooming) phenology of temperate and boreal trees (for more details see [1]). These models are based on the response of trees to the chilling requirement (the accumulation of which breaks endodormancy), and also the response of bud growth to forcing heat requirement (which forces the growth during spring after endodormancy break). Phenological models are important tools for planning agricultural practices [3], particularly for the prevention of frost damage in the short term [4,5] and for projecting the impact of climate change in the long term [6][7][8][9].
Dormancy is an important phase in the life of temperate zone fruit trees. Regular dormancy is essential for surviving unfavorable temperatures in winter [10]. Sufficient chilling units are needed for the breaking of endodormancy of apricot (Prunus armeniaca L.) buds. According to Andreini et al. [11], early cultivars are controlled by minimum temperatures, but intermediate and late cultivars are controlled by mean temperatures.
cultivars. In addition, we used the time when the string stage occurred in at least 50% of the trees [34]. Our collection was in Szigetcsép (47 • 15 10.6 N 18 • 57 26.6 E) until 2007. Three years before the excision of this plantation, we set up a new experimental plantation in 2004 in Soroksár (47 • 23 51.8 N 19 • 09 03.2 E), where we continued our studies. The environmental conditions of the two production sites are very similar, because they are very close to each other. The examination of the phenological process of apricot cultivars was carried out in parallel at both plantations in 2004, 2005, and 2006. No significant difference was found between the two production sites in these processes [34]. The beginning and end of the blooming of cultivars occurred on the same days at both sites. Data evaluated in this paper are from Szigetcsép between September 1994 and April 2006, and from Soroksár from September 2006 to April 2020. The training system at both plantations is compact vase with 5 × 3 m spacing with Prunus cerasifera seedling rootstocks. In the orchard, integrated plant protection was applied. In both plantations, six trees per cultivars were available for testing.
For the model building, we used the daily mean temperatures observed in the synoptic station of Marczell György main observatory of the Hungarian Meteorological Service (47 • 25 49.2 N 19 • 06 43.7 E) between September 1994 and April 2020. These data were considered to be good estimations of the daily mean temperature in the apricot plantations, taking into account the close distance between the observatory and experimental plantations and the low variance of the temperature layer, as well as the flat terrain of the surrounding area [54].

The Unified Model
The unified model [43] uses nine parameters (a c , b c , c c , b f , c f , w, k, C crit , t c ) to estimate the day of bud burst or the beginning of blooming of woody trees. In what follows, during the description of the unified model (UM), we use the definition of UM in the sense that it estimates the beginning of blooming.
The estimated critical amount of accumulated chilling units (C crit ) necessary to break the endodormancy and the total accumulated chilling units (C tot ) are defined as follows: where a c and b c are estimated weight parameters, c c is the estimated parameter of the maximum point, T is the daily mean temperature, t 0 is a fix date (1st of September), t 1 is the beginning of ecodormancy (beginning of the forcing period), and the estimated parameter t c is the last day of chilling unit accumulation. By fixing t 0 to 1st of September, we begin our chilling unit accumulation calculations on this date [43], but this does not constrain the model accuracy, because under Hungarian environmental conditions, chilling unit accumulation starts much later, mainly around the end of October [9,55].
The function of the chilling unit accumulation (1) or (2) is a bell-shape curve (Figure 1a) with its maximum value taken at the temperature that is optimal for the plant for chilling unit accumulation (c c ). This optimal value is usually a positive temperature close to zero. The estimated critical amount of forcing units ( ) necessary to stimulate the bud burst is calculated as follows: where is an estimated weight parameter, is an estimated parameter of the inflexion point of the function, is the daily mean temperature, is the beginning of ecodormancy (beginning of the forcing period), is the day of the beginning of blooming, > 0 is an estimated weight parameter, and < 0 is an estimated parameter of the relation of Fcrit and the total amount of chilling units Ctot. The beginning of blooming is set off when Ftot ≥ Fcrit. Considering Equation (4), it is obvious that Fcrit and Ctot are in a negative relationship [56]. This relationship expresses that a higher amount of chilling units means a lower amount of forcing units is needed to set the beginning of blooming off, which was proven experimentally [57][58][59][60][61]. The estimated critical amount of forcing units (F crit ) necessary to stimulate the bud burst is calculated as follows: where b f is an estimated weight parameter, c f is an estimated parameter of the inflexion point of the function, T is the daily mean temperature, t 1 is the beginning of ecodormancy (beginning of the forcing period), t 2 is the day of the beginning of blooming, w > 0 is an estimated weight parameter, and k < 0 is an estimated parameter of the relation of F crit and the total amount of chilling units C tot . The beginning of blooming is set off when F tot ≥ F crit . Considering Equation (4), it is obvious that F crit and C tot are in a negative relationship [56]. This relationship expresses that a higher amount of chilling units means a lower amount of forcing units is needed to set the beginning of blooming off, which was proven experimentally [57][58][59][60][61]. The function of the forcing unit accumulation (3) is an S-shape curve (Figure 1b), expressing the higher and higher forcing effect of warming temperatures.
We applied two simplifications to Chuine's unified model, because the high number of parameters can complicate the optimization, often impairing the accuracy of the estimate [62][63][64]. Moreover, there is a high probability that the underlying biological content disappears during optimization [7,65]. The simplifications are as follows: (1) According to Caffarra and Eccel [7], b c can be set to 0, so we described the c accumulation in the endodormancy with the following equation: (2) It does not strongly constrain the model accuracy if we assume that the chilling unit accumulation ends at the beginning of ecodormancy (t c = t 1 ), [43]. We defined this day as the one when the string stage occurs [66,67]. The observed string stage data were available from the data base of the examined apricot cultivars in the time period 1994-2020. Assuming t c = t 1 , it follows that C crit = C tot . In our study, t 2 is the day of the beginning of blooming.
We had six parameters left to estimate (a c , c c , b f , c f , w, k). Based on the literature, we chose a particularly large parameter space for the first runs [7][8][9]43,63,64]. The limits of this parameter space are shown in Table 1. Since k is a negative number very close to zero, we performed the optimization for k exp according to the following equation: Table 1. The boundaries of the parameter space and the step length 1 used for the first estimation.

Minimum Value Maximum Value
Step Length 0.01 1 These are not fixed step lengths: the σ parameter of the Gaussian distribution was used to select the step length.

Parameter Estimation with the Simulated Annealing Method
We used the simulated annealing method to estimate the parameters [51,52]. During the simulated annealing calculation process, we searched for the optimal six-dimensional parameter vector, which results in the best-fit to the measured blooming data set. The goodnessof-fit was defined as the root-mean-square error (RMSE) between the estimated and observed data. This measure was minimized within the boundaries of the parameter space. To achieve this, we started a random walk process within the chosen allowed parameter space ( Table 1). The starting coordinates were randomized according to a uniform distribution defined between the boundaries of the individual parameters.
At the beginning of the random walk, we set a so-called 'temperature' parameter T * to its starting value T start , defined as the goodness-of-fit at the starting point. In each random walking step, the temperature parameter T * was lowered ('cooled') by a factor of 1 − T step , where T step = 0.001. The 'cooling' process was performed from T start to T crit = 0.01·T start , so the random walk was finished when T * ≤ T crit .
At the beginning of each walk, we randomly chose a parameter to modify (i.e., the axis along which to step with optimization by 'cooling'), and a step length from a Gaussian distribution with zero mean and σ deviation (σ was fixed, Table 1). If a step moved out from the interval of allowed parameters, we stepped towards the opposite direction along the axis of the same parameter. Note that in some exceptional cases when this opposite step also moved out from the boundaries (i.e., the step length was too large), we terminated the whole fitting process, and sent an error message.
A step is accepted with a probability 1 if the measure of goodness-of-fit is smaller than it was previously, and it is accepted by a probability of p = e − Di f f T * if the measure is larger than before.
Here, Di f f is the absolute difference between the old and the new goodness-of-fit measure. One can see that as T * decreases in every step, the same Di f f results in a decreasing probability of the new position being accepted. This step acceptance algorithm ensured that the optimization is not stuck in local minima [53]. As a result of one optimization run (a 'walk'), we obtain the best-fitting vector of six parameters and the corresponding root-mean-square error (RMSE) value. Altogether, we performed 10,000 simulated annealing walks per cultivar to estimate the parameter vectors of the three apricot cultivars.
We used the MATLAB program system, Microsoft Excel (version: 2015), and Microsoft paint.net 4.2.15 software for our calculations and representations of our results.

Results
In reporting our results of 10,000 random walks to estimate the parameters of the three cultivars, we use the following notations for the cultivars: 'Ceglédi bíborkajszi' as 'cb', 'Gönci magyar kajszi' as 'gm', and 'Rózsakajszi' as 'ro '.
In what follows, we lead the reader through the optimizing procedure in eight steps.

1.
First, we observe that for each apricot cultivar, the forcing parameters c f and w are strongly related, although not linearly (see the upper panel in Figure 2). Therefore, these two parameters cannot be optimized independently. Using the empirical relationship between c f and w obtained from the optimization process, we calculate the optimal values of c f for all fixed values of w; Diversity 2022, 14, x FOR PEER REVIEW 6 of 15 the axis of the same parameter. Note that in some exceptional cases when this opposite step also moved out from the boundaries (i.e., the step length was too large), we terminated the whole fitting process, and sent an error message. A step is accepted with a probability 1 if the measure of goodness-of-fit is smaller than it was previously, and it is accepted by a probability of = * if the measure is larger than before. Here, is the absolute difference between the old and the new goodness-of-fit measure. One can see that as * decreases in every step, the same results in a decreasing probability of the new position being accepted. This step acceptance algorithm ensured that the optimization is not stuck in local minima [53]. As a result of one optimization run (a 'walk'), we obtain the best-fitting vector of six parameters and the corresponding root-mean-square error (RMSE) value. Altogether, we performed 10,000 simulated annealing walks per cultivar to estimate the parameter vectors of the three apricot cultivars.
We used the MATLAB program system, Microsoft Excel (version: 2015), and Microsoft paint.net 4.2.15 software for our calculations and representations of our results.

Results
In reporting our results of 10,000 random walks to estimate the parameters of the three cultivars, we use the following notations for the cultivars: 'Ceglédi bíborkajszi' as 'cb', 'Gönci magyar kajszi' as 'gm', and 'Rózsakajszi' as 'ro '.
In what follows, we lead the reader through the optimizing procedure in eight steps.
1. First, we observe that for each apricot cultivar, the forcing parameters and are strongly related, although not linearly (see the upper panel in Figure 2). Therefore, these two parameters cannot be optimized independently. Using the empirical relationship between and obtained from the optimization process, we calculate the optimal values of for all fixed values of ;

2.
Based on the resulted parameter vectors of all walks, we plot the histogram of the optimal parameter values of w (see the lower panel in Figure 2). We see that the optimal parameter values of w are dense mainly around three or four small to high values; more exactly, around one small, two medium, and one high value for 'Rózsakajszi C.1406', while around one small, one medium, and one high value for the other two cultivars. We immediately exclude the high values, because we obtained C f values in those cases between −10 • C and −30 • C, which are unlikely during the forcing period in Hungary [9,68]; 3.
The results of most walks are dense around the small value for each apricot cultivar, and we obtained the lowest RMSE values here, too. So, we fix the parameter value w at the median of the preferred range of 'small' optimal parameters w: w cb = 15.6, w gm = 17.9, w ro = 19.4; 4.
As a next step, we narrow the parameter space according to the biologically possible parameter values for Hungary (Table 2) [7,9,65,69,70]. In the original parameter space, we find several similarly good parameter vectors, that fit statistically very well to the observed blooming dates, but they are biologically impossible. Table 2. The boundaries of the parameter space and the step length1 for the second estimation.

Minimum Value Maximum Value
Step Length With the aim of finding biologically realistic optimal parameter vectors [68], we started another 10,000 walks for each cultivar, applying the fixed parameters w cb , w gm , and w ro . As a result, we obtained 10,000 limit vectors of the six parameters that the optimization converged to. In order to visualize this result, we created fifteen histograms (for five parameters of the three cultivars based on 10,000 limit values, not shown) to find the local optimum parameter values. To detect the significantly more frequent limit values, we divided the histograms into one hundred bins, because it is the square root of the amount of the limit parameter values. Thus, in the random case, an average of one hundred data points falls in each bin (µ = 100). Unsurprisingly, this happened with natural fluctuation following a Poisson distribution. Having such a large number of limit points, the Poisson distribution could be approximated by the Gaussian distribution with standard deviation σ = √ µ = 10. Focusing on a histogram, we define the limit point frequency in a bin as significantly high if it is above the value µ + 3σ = 130. From the optimized 15 parameters, in seven cases, namely for c cb f , k cb exp , a gm c , c gm c , c gm f , a ro c , and c ro f , we obtain one or more significantly frequented bins next to each other. In another six cases, namely for the parameters b cb f , b gm f , k gm exp , c ro c , b ro f , and k ro exp , there are more non-adjacent significantly frequented bins. In the remaining two cases (a cb c , c cb c ), we do not obtain any bin with at least 130 data points, but there are several bins with more than µ + 2σ = 120 optimized parameter values. We define the local optima of the parameters as the medians of the significantly frequented bins. In the case of non-adjacent significant bins, we searched for the local optimum of the parameter in the bin with the largest number of optimized values (Table 3). Table 3. Characterization of the local optimum parameter vectors of the unified model: the minimum and maximum values of the bins with significantly higher number of optimized limit values (optimum bins), the number of the limits in them, the median and the standard deviation, together with the maximum, the median, and standard deviation of the RMSE. The values belonging to the cultivars 'Ceglédi bíborkajszi', 'Gönci magyar kajszi', and 'Rózsakajszi C.1406' are denoted by the superscripts 'cb', 'gm', and 'ro' respectively.  , and c ro f ), the global optimal parameter values do not fall in the local optimum bins. This is most surprising for the parameter c ro f , where more than 70% of the walk limits fall in the local optimal bin, but the global optimum parameter value does not;

6.
Using the global optimum parameter vector, we estimated the blooming date for each apricot cultivar with an average error less than 2.5 days (RMSE < 2.5). For comparison, if we take the mean blooming data calculated over all the years as a constant [64], the average error of the estimation (i.e., the error of the base model) is as high as 9.7-10.6 days, depending on cultivars; 7.
Finally, based on the daily average temperature and the global optimal parameter vectors, we calculated the critical amount of chilling and forcing units, and determined the chilling and forcing process for each cultivar in the period 1994-2020 (Figures 1 and 3). We provide the parameter values that are optimized with the simulated annealing method, applying the unified model and the observed string stage and blooming data of years 1994-2000 (Table 4). The temperature that is optimal for the plant for chilling unit accumulation (c c ) is 1.50 • C for 'Ceglédi bíborkajszi', 2.13 • C for 'Gönci magyar kajszi', and 2.42 • C for 'Rózsakajszi C.1406' in the period 1994-2020. According to our calculations, the most chilling units (C crit = 29.8 units) are necessary for 'Gönci magyar kajszi', and the least chilling units (C crit = 12.7 units) are required by 'Ceglédi bíborkajszi' for breaking the endodormancy ( Table 4). The inflection point of the forcing unit accumulation (i.e., c f ) is between 8.30 and 9.04 • C, depending on the cultivars. This curve has no maximum point, but the forcing unit accumulation is close to the maximum (1 unit) at 12-15 • C (more than 0.9 units) that could be considered as 'optimal temperature' for the plant in their preparation for blooming. The average accumulated forcing units for the blooming are between 14.0 and 16.4 units for each cultivar in the period 1994-2020 (Table 4). Surprisingly, the absolute value of parameter k of our results is larger than is reported in the publications of other researchers (i.e., in between −10 −4 and −10 −8 ) [8,9,43,64]. This may lead to a conclusion that, in the case of Hungarian apricots, the chilling unit accumulation has a relatively larger effect on forcing unit accumulation. Note that though the change in the observed and estimated blooming start is proven as significant in the past 26 years, the slopes of their trends are all negative, bel −0.2 (Table 5). Since the slopes of the parameters , = , and are alm equal to zero, i.e., they seem to not be changing with time as the blooming start does, model can be directly applied in climate change impact studies as well. While the un isfactory or slowed down chilling unit accumulation delays the blooming start, the wa ing temperatures speed up the forcing accumulation. This complex relation with oppo directions can be simulated by our model. Note that though the change in the observed and estimated blooming start is not proven as significant in the past 26 years, the slopes of their trends are all negative, below −0.2 (Table 5). Since the slopes of the parameters F crit , C crit = C tot , and F tot are almost equal to zero, i.e., they seem to not be changing with time as the blooming start does, our model can be directly applied in climate change impact studies as well. While the unsatisfactory or slowed down chilling unit accumulation delays the blooming start, the warming temperatures speed up the forcing accumulation. This complex relation with opposite directions can be simulated by our model.

Discussion and Conclusions
Phenology models play an important role in planning agricultural practice [3][4][5], and also in research of phenological changes in plants caused by climate change [6][7][8][9]. In our work, we fit Chuine's unified model [43] to the blooming data sets of three apricot cultivars ('Ceglédi bíborkajszi', 'Gönci magyar kajszi', and 'Rózsakajszi C.1406') in Hungary in time period 1994-2020. Szalay et al. [34] analyzed the same data sets in the period 1994-2018. They observe a shift in string stage of ca. 0.5 days/year later, and in the beginning of blooming of ca. 0.125 days/year earlier. They explained this change by the warming climate in the recent decades in Hungary [71]. The chilling units accumulate more slowly in mild winters, which delays the phenological phases in dormancy [72,73]. On the other hand, the warmer spring causes an earlier beginning of blooming [31].
Due to the higher winter temperatures, and the inhibitory effect on chilling unit accumulation, the phenology models that consider only forcing unit accumulation are no longer suitable for estimating the beginning of bud burst and blooming [74][75][76][77]. That is why we chose the unified model [43], which contains both chilling and forcing unit accumulation. For the estimation of the parameters, we used the simulated annealing method [51,52], which avoids the problem of getting stuck in local minima.
We reduced the nine-parameter model to a six-parameter model with two simplifications [7,43], thus, we could accelerate the parameter estimation [63,64]. We show that there is a strong, non-linear relationship between the parameters c f and w. After obtaining a great number of parameters from our random walk optimization, using the histograms of the limit parameters, we defined the local optimum parameter values as the median of the bins with a significantly high number of limits that the optimization converged to. Finally, we defined the global optimum parameter vector as the one with the lowest root-mean-square error (RMSE). It is an interesting result that two-thirds of the global optimum parameter values do not fall into the local optimum parameter bins.
We estimated the beginning of blooming during the period of 1994-2020, using the global optimum parameter vector with an RMSE of less than 2.5 days, which can be considered as a good estimation compared to the results of several papers (2.6-5.6; [8,9,68,78]).
In our calculations, the maximum of the chilling unit accumulation is between 1.50 and 2.42 • C, depending on the cultivars. The inflection point of the forcing unit accumulation curve is between 8.30 and 9.04 • C, while the forcing unit accumulation is almost maximal (more than 0.9 units) between 12 • C and 15 • C, depending on the cultivars. Other research [68,79,80] determined the maximum of the chilling and forcing unit accumulation for European fruit trees at temperatures that are slightly higher: around 5 • C and around 25 • C, respectively. But there may be differences in response to temperature even within populations of a given species [81]. Thus, it is possible that the results of various studies cover a relatively wide range of chilling and forcing unit requirement values. The maximum of chilling unit accumulation is found between −28 • C and 9 • C for Quercus mongolica [63] and Taxus baccata [43], respectively, and the inflection point of forcing unit accumulation is shown to be between −14.5 • C and 16 • C in the case of Olea europaca [43] and Vitis vinifera L. [7], respectively. In the case of extreme negative values, the limit is probably statistically the best parameter vector, but biologically it may be incorrect or hardly interpretable [68]. In the case of mechanistic phenology models, however, the biological representativeness is essential.
In some cases (when parameter k is around −10 −7 or more), a simplification can be applied with k = 0 and w = F crit . But, in view of our research, this simplification is not justified, because the absolute value of parameter k is larger (between −0.0072 and −0.0086) than it is in publications of other researchers (−10 −8 -−10 −4 ) [8,9,43,64]. It can mean that, in the case of Hungarian apricots, the chilling unit accumulation has a larger effect on forcing unit accumulation.
With our adapted Chuine's unified model, in the future we can provide a more accurate estimation for the blooming time of the apricot varieties in the studied area. Therefore, the preparation for frost-hazardous conditions can become notably more effective. Since the blooming times of apricot varieties in Hungary are quite close to each other [21], and the climatic conditions of apricot growing areas are also very similar [82], our results can easily be adapted to other Hungarian varieties and regions.
Making estimations for the future using climate model data requires other preliminary studies, because the extension of the phenological models for the future may be questionable [7,9,64,68,83,84], though we also can find very promising new results [85,86].