Gamma, Gaussian and logistic distribution models for airborne pollen grains and fungal spore season dynamics

The characteristics of a pollen season, such as timing and magnitude, depend on a number of factors such as the biology of the plant and environmental conditions. The main aim of this study was to develop mathematical models that explain dynamics in atmospheric concentrations of pollen and fungal spores recorded in Rzeszów (SE Poland) in 2000–2002. Plant taxa with different characteristics in the timing, duration and curve of their pollen seasons, as well as several fungal taxa were selected for this analysis. Gaussian, gamma and logistic distribution models were examined, and their effectiveness in describing the occurrence of airborne pollen and fungal spores was compared. The Gaussian and differential logistic models were very good at describing pollen seasons with just one peak. These are typically for pollen types with just one dominant species in the flora and when the weather, in particular temperature, is stable during the pollination period. Based on s parameter of the Gaussian function, the dates of the main pollen season can be defined. In spite of the fact that seasonal curves are often characterised by positive skewness, the model based on the gamma distribution proved not to be very effective.


Introduction
The characteristics of a pollen season, such as timing, curve and magnitude, depend on a number of factors such as the biology of the plant, species spectrum and environmental conditions. Patterns in the occurrence of airborne pollen are largely chaotic, but some features of pollen seasons of specific taxa are characterised by relatively low variability (Wołek and Myszkowska 2008). A similar opinion was presented by Hyde (1956) who stressed that for many pollen taxa, it is possible to determine a ''dense'' period during which the pollen is present in the air at its maximum concentrations every year, e.g. for Ulmus or Fagus. Time series of pollen seasons from the same taxon in successive years can follow a different pattern (Comtois and Sherknies 1991;Piotrowska and Kubik-Komar 2012). This phenomenon is more noticeable in taxa with lower airborne pollen concentrations such as, for example, Rumex. This simply results from the statistics of pollen counts.
Moving averages, with varying degrees of smoothing, can be used to represent the curve of the pollen season. In the case of many taxa, this method does not always produce the expected result and it is still difficult to discern any trends or the characteristics of the season. When describing Urticaceae pollen time series, Belmonte and Canela (2002) applied a nonparametric method, Friedman's super smoother. This method chooses from among three smoothers by cross-validation for the best prediction.
Mathematical models describing pollen season curves have also been presented in the literature. They are a greater or lesser simplification of the reality, but they provide a better understanding and interpretation of aerobiological phenomena and allow these phenomena to be predicted. The occurrence of airborne pollen/fungal spores has a seasonal nature, and the curves representing their temporal concentrations often have a shape that is close to a bell curve. A model based on a Gaussian function has been proposed by Kasprzyk (2011). On the basis of the function parameters, the author defined the date of maximum concentration and pollen season duration. This method is not only useful in the case of a unimodal curve of the pollen season, Kasprzyk and Walanus (2010) and Belmonte and Canela (2002) have fitted this type of model to multimodal pollen seasons.
Pollen grains and fungal spores can remain airborne for a long period of time after emission. The post-peak period is often longer than the pre-peak period, and the curve of the season is often characterised by positive skewness. As a result, gamma distribution-based models have been successfully used to describe the time sequence of pollen and fungal spores by Comtois (2000) and Belmonte and Canela (http://lap.uab.cat/ aerobiologia/general/pdf/altres/MCC_XIIPC_slides. pdf; http://lap.uab.cat/aerobiologia/general/pdf/altres/ TESAGamma.pdf), such models are successful in describing the seasonality of airborne pollen and fungal spores.
Another universal model has been presented by Ribeiro et al. (2007). It is a logistic model and its effectiveness does not depend on the duration or annual variation of the pollen season, or the taxa. By its nature, the model adapts itself to the cumulative curve of pollen concentration, as proposed by Pathirane (1975). The authors tested the model on several taxa, and it always produced an exceptionally good fit, almost 99 %. They stressed that the pattern followed by pollen seasons of various taxa from different regions can be quickly compared by the means of the function coefficients, in particular the growth rate coefficient. The models based on this function are often used to describe growth processes of living organisms (Gregorczyk 1991). A different approach was used by Stepalska and Wołek (2005). The period of occurrence of airborne Alternaria spores is exceptionally long and irregular. They divided it into two time intervals: the pre-peak and post-peak period, and both periods were fitted with a mathematical function that best described it. The best fitting was found for the exponential function, while the fitting for the power function was lower. Other models have been proposed that are based, among other things, on Poisson regression and Weibull distributions, gradient boosting method, transfer function model and autoregressive model (Limpert et al. 2008;Ocana-Peinado et al. 2008;Hilaire et al. 2012).
The main aim of this study was to develop mathematical models that explain dynamics in atmospheric concentrations of pollen and fungal spores. Plant taxa with different characteristics in the timing, duration and curve of their pollen seasons, as well as several fungal taxa, were selected for this analysis. Several models were examined, and their effectiveness in describing the occurrence of airborne pollen and fungal spores was compared. The effectiveness of the models was tested for on data collected at two separate sites differing in geobotanical and climatic conditions.

Study area
Aerobiological monitoring was carried out in the centre of Rzeszów (N50°01 0 45 00 , E22°00 0 57 00 ). Three years (2000)(2001)(2002) were arbitrarily selected from several years of the database. Rzeszów is located in south-east Poland in the province of Carpathian Foothills. Mean annual precipitation is above 730 mm, and mean annual temperature is 8.8°C. The warmest month is July (17.5°C), and the coldest is January (-4.6°C). The average vegetation period is 215-220 days. Rzeszów city boundaries have an agricultural character; the environs are a mosaic of forests and crop fields. The obtained results were compared with ones achieved in two other localities. The authors only possessed data from 1995.
• Zakopane (900 m a.s.l.) is a small town located in the Tatra Mountains, southern Poland. The mean annual precipitation ranges from 900 to 1,300 mm. The mean temperature in January is -4.9°C and in July 14.7°C. The average vegetation period is 180 days. In the close vicinity of Zakopane, there are more meadow areas than farmland. Fagus silvatica, Abies alba and Picea excelsa are an important component of forests. • Ostrowiec Ś więtokrzyski (175 m a.s.l.) is a midsized city in Małopolska Upland in central Poland, where the mean annual precipitation is 550-650 mm and the vegetative period is 210 days. The mean temperature in January is -3.5°C (the coldest month) and in July 18°C (the warmest month). Vegetation in the town and its surroundings consists of ruderal vegetation, pine forests, semi-natural community of grasses and anthropomorphic habitats.

Aerobiological monitoring
Investigations were continuously conducted at all sites using volumetric spore traps of the Hirst design (Hirst 1952). In Rzeszów, the trap was located about 12 m above the ground level, in Ostrowiec Ś w. 36 m and in Zakopane 10 m. Analyses by light microscope were carried out at magnification 4009. Pollen grains were counted along 12 vertical transects on each microscope slide. Each transect corresponds to a 2-h interval. The obtained result was expressed as the number of pollen grains in a cubic metre of air daily average (grains/m 3 ) (Frenguelli 2003). Three pollen taxa with different characteristics in the timing, duration and curve of their pollen seasons were selected for analysis: Artemisia, Betula, Rumex. In addition, four fungal taxa were taken into consideration: Cladosporium, Alternaria, Ganoderma and Botrytis.

Gaussian model
The bell curve fitting method [Gaussian curves, exp (-x 2 )] was used to describe the seasons of airborne pollen and fungal spores.
where y, daily count of pollen grains/fungal spores m 3 / 24 h; t, day from the 1st January; t 1 , day of the modelled maximum pollen/spore concentration expressed as the number of days from the 1st January; a, the value of the modelled maximum pollen grains/ spore concentration m 3 /24 h; and s, SD (in the sense of normal distribution) as a measure of half of duration of the maximum pollen/fungal spore season. Using this function, a new method for defining the seasonal occurrence of airborne pollen and fungal spores was proposed. In the case of the Gaussian curve, 95.4 % of all values are theoretically included within the range of the t 1 ± 2s. The dates calculated on this basis were compared with the dates based on the 95 % method (Goldberg et al. 1988), whereby the start of the season is the day when the cumulative sum of airborne pollen and fungal spores is 2.5 % of the total sum, and the end of the season is when this value reaches 97.5 %. This method for defining the limits of the season is frequently used in aerobiological studies (e.g. Stach et al. 2008).
Percentiles were also determined using the Excel formula for inverse normal distribution. Using the parameters determined on the basis of fitting this function, the occurrence of airborne pollen and fungal spores was presented as a cumulative concentration.

Gamma model
Gamma distribution that includes the right skewness of the curves representing the occurrence of pollen and fungal spores in the atmosphere.
a, the value of the maximum pollen grain/fungal spore concentration per m 3 /24 h; c, shape parameter (the lower c, the larger skewness); s, scale parameter; t, number of days from the 1st January; and t g , shifting parameter for t 1 .
Day of the maximum pollen/fungal spore concentration can be expressed as a function of the above parameters: On the basis of this model, the dates of seasons of airborne pollen and fungal spores were defined using the inverse gamma distribution function in Excel.

Logistic differential model
The logistic function also has its differential version, which is less known. The formula is given below. The curve of the function can be normalised to the value of 1 as its maximal value (the same is necessary in the case of Gaussian distribution). In addition, if the horizontal axis (time axis) of the function will be normalised to obtain the same time space, e.g. 68 % quantiles, the two differential functions (the Gaussian and the logistic) produce very similar results. There are two different functions: the main difference between them is the way that the functions approaching zero. The Gaussian function reaches zero quickly, like the function exp(-x 2 ), while the logistic function is simply the exponential in that aspect; exp(-x). Whether this difference is essential for modelling pollen data is a question that needs to be answered experimentally. Comparing the results of fitting the integral and differential functions (to the integral and differential data) is not straightforward. The formula of logistic differential function is: Þ=s t, number of days from the 1st January; t 1 , day of the modelled maximum pollen/ fungal spore concentration expressed as the number of days from the 1st January; a, the value of the modelled maximum pollen grains/fungal spores concentration m 3 /24 h; and s, SD as a measure of half of duration of the maximum pollen/fungal spore season.

Cumulative models: Gaussian and logistic
Time series of the pollen grains/fungal spores collected can be visualised/treated in two different ways. The first treatment is strictly connected to the method of data collection. The number of pollen grains/fungal spores is counted every day. That natural way can be called differential, since we do not simply have a number of items (pollen grains), but the number of items calculated (counted) per unit of time (day). That starts from zero on the 1st January and approaches a maximum somewhere in the summer (e.g.) and reverts to zero in December. Such behaviour in time is explained by, e.g. the bell-shaped curve (Gaussian probability distribution).
Another cumulative function frequently used for a wide variety of different purposes is a logistic one.
a, total sum of airborne pollen and fungal spores; t, number of days from the 1st January; b, position parameter; and c, scale parameter.

Rejection of outliers
Pollen and fungal spore data, like data of any origin, are not considered to be especially resistant to outliers.
In the example (Fig. 1a), two outliers are evident in the post-peak period of the season. There is the possibility of simple gross errors, mistakes created when producing and processing the data. The most evident outlier in Fig. 1a is surrounded by days with low, but nonzero counts. Such points may be treated as outliers as well, depending on the given taxon and its behaviour of pollen season. The models (not only Gaussian) can be considered robust where outliers far from the main season are concerned. This is because the function is close to zero in such a region. As a result, the sum of squares of differences is weekly dependent on small changes in the model parameters. A single outlier is too weak to move the whole model against the bulk of the 'healthy' data. As such, the simple idea of least squares works well with pollen data.
Nevertheless, evident outliers should be removed. We propose the simple rule, based on the 'basement' of the four surrounding days. The outlier (x i ) should be rejected and replaced with the average value of neighbours when where sum is a total, yearly sum of pollen grains/ fungal spores. The value of 5 in the formulae is chosen in an arbitrary way. The number 100 in the first part of the conjunction is arbitrary as well. However, it is of less influence. That part of the conjunction removes small counts from the suspicious range, what is necessary if the second, the main one part, of the conjunction is based on the relation to the surrounding days that may be empty of pollen (Fig. 2). One pollen grain surrounded by no pollen days cannot be treated as an outlier (1 [ 590).
The so-called Friedman's 'super smoother' was used to visualise time series of airborne pollen and fungal spores. It has been applied with statistical computing environment R (www.r-project.org). It is a 'running lines' smoother that chooses between three spans for the lines based on minimising cross-validated error.
The function 'nonlinear estimation' Statistica 9.0 software was used for the analysis. The Levenberg-Marquardt method was applied to determine the initial values. The SE was given for all the function parameters. The goodness of fitting of the presented functions was expressed by the coefficient of determination. The statistical significance of the parameters of all functions was accepted at a B 0.05.

Models for pollen grains
In the case of the pollen taxa analysed, all the Gaussian function coefficients were statistically significant, but the goodness of fit of the function to the actual values varied significantly. The goodness of fit was also dependent on the type of taxa. The multi-species Rumex with medium or low concentrations of airborne pollen is characterised by poorly fitted functions. The Gaussian function was able to describe the occurrence of Artemisia pollen in the air in both Rzeszów and Ostrowiec Św. (Fig. 1a). The coefficient denoting the estimated maximum value showed the largest error. The other coefficients were characterised by lower SEs (Table 1). The function based on the Gaussian distribution appears to be resistant to outliers. This can be seen on the example of Artemisia in 2000. In that year, two jumps in concentrations were recorded after the main period of pollen release (Fig. 1a). These outliers would generally be expected that these outlying points would have some influence on the model, but this did not occur. Replacing outliers by the average of the surrounding days changes the values of the model and lowers the SE of estimation ( Table 2).
The gamma function-based model shows the poorest fitting for pollen seasons. For many taxa, it was not possible to develop a model whose coefficients would be statistically significant ( Table 6). As seen in Fig. 1a, the curve for Artemisia in 2000 only fits well in the post-peak period (the fitting for the pre-peak period was very poor). This is characteristic of all the taxa analysed. As in the case of the Gaussian function, a new model was prepared for Artemisia in 2000 after outliers were rejected. The goodness of fit of the model to the data was slightly better (it increased from 76 to 83 %), but the SEs for the function parameters were higher (Table 3).
The fitting of the differentiable logistic function produced good results, comparable with the Gaussian function (Tables 1, 4). In accordance with the principle adopted by the authors for Artemisia 2000, outliers were rejected and this model was refitted. The goodness of fit increased from 76 to 83 % (Table 4).
The logistic function fitted the cumulative curves of pollen seasons of all the taxa analysed almost perfectly (Table 5; Fig. 1b). For Artemisia, the values of the parameter c were the highest and similar in successive years. This indicates a rapid increase in airborne pollen concentrations during the season. In the case of Rumex, the values of this parameter were the lowest, which indicates an entirely different pattern of the pollen season. There was a slow rate of increase in pollen grains of this taxon in each year of the study and at each site (Table 5).

Models for fungal spores
It was possible to fit the Gaussian distribution-based models to each type of spore in each year for all sites. Plotted are days only which confirm the criterion of share more than 1 % of the total sum. The plotted value (vertical axis) is that used in the criterion: Þ . The value 5 seems to be good limiting one. The case of Artemisia concentrations in 2000 in Rzeszów is presented The goodness of fit ranged from 23 % for Botrytis in 2002 to 82 % for Cladosporium in 2001. Although all the coefficients of the model were statistically significant, the highest error related to the coefficient a, which denotes the estimated maximum value. The errors of the coefficients relating to the peak date (t 1 ) and to the parameter s were much lower (Table 1). It was not possible to fit the gamma distributionbased model in all cases. In a few cases, the goodness of fit of the gamma distribution-based function was a slight improvement compared with the model based on the normal distribution. This function was the best at describing the seasonality of occurrence of airborne Botrytis spores in 2001, while the poorest fit was found also for Botrytis in the following year. In the case of Cladosporium, the curve of the gamma model exhibits a very good fit for only the post-peak period (Fig. 3a). The shape parameters indicating the skewness of the time series had similar results. For example, high values ([3) showed that the Cladosporium in Ostrowiec Świętokrzyski seasons had low skewness. Low values of these coefficients, e.g. Cladosporium in Zakopane and the Ganoderma season in 2000 in Rzeszów, indicate that these spores remained in the air for a long period of time at low concentrations ( Table 3).
The models based on a differentiable logistic function showed a rather good fitting in many cases. The fitting was only below 50 % in four cases. The majority of R 2 coefficients were slightly higher than those for the Gaussian function, while the errors of the parameters a and t 1 were slightly lower. This model fitted much better to the data in the case of Ganoderma in 2000, Botrytis in 2001 (Table 4).
The logistic function produced an almost perfect fit for the rate of increase in the number of airborne fungal spores in all cases (Fig. 3b). However, parameter a, which denotes the cumulative total of spores, had a large or very large error. The values of parameter c are generally very low, which shows that the rate of increase in concentrations was also very low ( Table 5). The logistic model curve roughly fitted the cumulative curve for Cladosporium spores in 2001 (Fig. 3b).

Graphical presentation of seasons of airborne pollen and fungal spores
All presented methods smooth the diagrams presenting the occurrence of airborne pollen or fungal spores to a greater or lesser degree. After visual analysis of Figs. 1a, 3a, we state that the Friedman's 'super smoother' method is very good. It described the unimodal season as well as multimodal seasons like Cladosporium in 2001.

Calculation of dates of pollen/fungal spore seasons
The use of a Gaussian function smooths the curves in airborne pollen and fungal spore concentrations. Sporadic occurrences of airborne pollen and fungal spores, i.e. days when pollen or fungal spores were recorded before or after the period of continuous occurrence, were automatically rejected from the model. For example, there were two episodes recorded in 2000 when airborne Artemisia pollen increased for just 1 day; at the end of August and at the beginning of October. Pollen seasons defined using the Gaussian and gamma functions failed to capture these peaks. The models determined the end of the season to be: (1) Gaussian function-end of August; (2) gamma function-beginning of September. According to the traditional percentage method, the season ended at the beginning of October (when 97.5 % had been reached) and both peaks were included. The two methods for defining pollen seasons based on the Gaussian function were then compared. The pollen season defined using percentiles from the modelled data was always shorter than the one defined using the method t 1 ± 2s. The differences in end date were either very small or there were no differences, for instance for Artemisia in 2000 and 2002 as well as for Alternaria in 2002. Much larger differences were found in the case of start date-the dates based on percentiles were always earlier than for the method t 1 ± 2s (Table 6). The relative SE of estimation of the model parameters (SE est ) is also given, for comparison. Since the SE out is much less than SE est the removing of outliers or not is simply irrelevant. Model is the Gauss function (a 1 , t 1 , s parameters of Gauss function)

Discussion
Aerobiological data are an example of data that do not have a normal distribution (Comtois 2000). They are characterised by positive skewness, and different data transformations (finding the logarithm, extraction of roots) do not change their basic shape. This limits the range of methods that scientists can use to analyse these data (Stępalska and Wołek 2005;Fernandez-Llamares et al. 2013). However, if we illustrate the occurrence of airborne pollen and fungal spores in time, it can be seen that the curves are close to being symmetrical and it is even clearer after Fisher's smoothing method is applied (Belmonte and Canela 2002;Hilaire et al. 2012). A question therefore arises whether a good model for defining seasons of occurrence of airborne pollen and fungal spores can be prepared based on the Gaussian distribution and the differentiable logistic function (unimodal functions). According to the authors' knowledge, the latter function has not previously been used for this type of data. The Gaussian and differential logistic models shown in this paper were very good at describing pollen seasons with just one peak. These are typical for pollen types with just one dominant species in the flora and when the weather, in particular temperature, is stable during the pollination period. Such types of season are characteristic of Artemisia. According to a floristic study conducted in Rzeszow by Ś więs (1993), the most frequently found mugwort species is A. vulgaris. Modified Gaussian and differential logistic models can produce a good fit for multimodal pollen seasons in which concentrations are relatively high and the peaks are clearly defined. Such a model for the Gaussian function was presented by Kasprzyk and Walanus (2010) for a bimodal pattern of the grass pollen season where it performed well (on average with 70 % explanation). Belmonte and Canela (2002) successfully used this model for the multimodal Urticaceae season. The Gaussian and differential logistic functions poorly describe pollen seasons of taxa, which have: pollen that remains in the air for a long time; low levels of pollen in the air; and peaks that are not particularly evident. Such seasons are characteristic of the multi-species genus Rumex (Piotrowska 2012). The functions do not fit the data particularly well, generally at a level lower than 50 %. Such a result is not satisfactory. Alternaria and Cladosporium spore Table 3 Parameters of gamma functions for chosen taxa (only statistically significant) 191.0 ± 2.8 64 ± 17 0.86 ± 0.13 The formulae y concentrations are high or very high and the peaks are often clearly defined, but the Gaussian and differential models often failed to fit the data. The seasonal pattern is strongly affected by external factors, such as the availability of substrates and weather conditions, as well as agronomic treatments that could alter the amount of fungal spore in the air (Corden et al. 2003;Skjøth et al. 2012). The goodness of fit is not dependent on where research is carried out, but the behaviour of the pollen season and type of taxa. The curves representing the occurrence of pollen grains of some taxa are often characterised by a longer post-peak period (Wołek and Myszkowska 2008;Myszkowska et al. 2011;Dahl et al. 2013;Fernandez-Llamares et al. 2013). Grewling et al. (2012) report that although Betula pollen seasons can last more than 1 month, 85 % of birch pollen appears within the first 2 weeks, and if the temperature reaches 20°C, it is 90 % within the first few days. Myszkowska and Piotrowicz (2009) state that Betula pollen concentrations are lower and the pollen seasons are longer when there are unstable weather conditions during the period preceding pollen release; in such cases an increase in pollen concentration was often recorded at the end of the season. The analysis of our results showed a similar pattern of the seasons of some taxa. Therefore, an attempt was made to develop a model based on the gamma distribution, which should incorporate these features. Belmonte and Canela indicated that this was a good method both for pollen and fungal spore seasons (http://lap.uab.cat/ aerobiologia/general/pdf/altres/TESAGamma.pdf; http:// lap.uab.cat/aerobiologia/general/pdf/altres/MCC_ XIIPC_slides.pdf). The presented results do not confirm this. The models have a poor fit or they are not statistically significant. Despite the fact that season dates can be defined on the basis of the function parameters, it seems that they do not correspond to reality. The beginning of the season is defined too late. In the case of Artemisia, this date only differs from the peak date by 3 days.
Different methods presenting cumulative data are given in the literature. The cumulative version of the distribution function visualises the main feature of the time series, i.e. the maximum of the season, in a rather hidden way. Cumulative curves are S-shaped and such a curve can be produced for the growth of organisms or, as proposed by Pathirane (1975) and Ribeiro et al. (2007), for the pattern of airborne pollen concentrations. They can be described using models based on the logistic function, and its basic formula and numerous modifications have been presented, among others, by Gregorczyk (1991) for plant growth. The presented results confirm the statement of Ribeiro et al. (2007) that the model performs well and the strength of the fitting is not dependent on the type of taxon. Compared with logistic functions, the cumulative Gaussian distribution performed poorly. Outlying or suspicious values (a mistake in the identification of pollen grains, equipment error) can reduce the effectiveness of models. Decisions about the rejection of atypical or suspicious results should be made on a case by case basis, taking into account all possible factors affecting the pattern of occurrence of airborne pollen and fungal spores such as: taxon, the specificity of pollination or sporulation, the flora of the region, human activity, or the purpose of the research.
Several methods for defining the limits of the pollen season are given in the literature. Jato et al. (2006) compared them in detail and proposed that the type of taxon (mono-or multi-species) and sensitivity to weather changes should be taken into account in defining pollen season dates. We have defined the seasonal occurrence of airborne pollen and fungal spores by making an assumption that they are regular and close to symmetry. We have based the definitions on the Gaussian distribution what is novelty in aerobiological study. The dates of seasons can be very quickly calculated having only two parameters: the estimated time of maximum and SD (t 1 ± s). The proposed method should be used when we are interested in a condensed season. The second method is based on percentiles calculated from adjusted results, and it is more time-consuming. The analysis of the results shows that the type of method should be adjusted to the type of taxa, which is consistent with Jato et al. (2006) suggestions.

Conclusions
The models based on the Gaussian and differentiated logistic functions are robust and suitable for pollen grains and fungal spores. They are recommended for unimodal seasons. They are quite resistant to outliers. Two novel methods for defining pollen seasons have been proposed for aerobiological studies. They are based on the Gaussian function. We suggest that they should be used to define the main pollen season. Comparisons of season dates defined by several methods show once again that, depending on the purpose of the research, the use of the same method for all taxa is not always justified and; moreover, one should consider whether in some cases the start and end of the season should be defined using different methods. In spite of the fact that seasonal curves are often characterised by positive skewness, the model based on the gamma distribution was not very effective. Friedman's smoother method is very good for graphically presentation of time series of pollen and fungal spores. It is not easy method; it needs the knowledge of R programme language.