Machine Learning for Predicting Mycotoxin Occurrence in Maize

Meteorological conditions are the main driving variables for mycotoxin-producing fungi and the resulting contamination in maize grain, but the cropping system used can mitigate this weather impact considerably. Several researchers have investigated cropping operations’ role in mycotoxin contamination, but these findings were inconclusive, precluding their use in predictive modeling. In this study a machine learning (ML) approach was considered, which included weather-based mechanistic model predictions for AFLA-maize and FER-maize [predicting aflatoxin B1 (AFB1) and fumonisins (FBs), respectively], and cropping system factors as the input variables. The occurrence of AFB1 and FBs in maize fields was recorded, and their corresponding cropping system data collected, over the years 2005–2018 in northern Italy. Two deep neural network (DNN) models were trained to predict, at harvest, which maize fields were contaminated beyond the legal limit with AFB1 and FBs. Both models reached an accuracy >75% demonstrating the ML approach added value with respect to classical statistical approaches (i.e., simple or multiple linear regression models). The improved predictive performance compared with that obtained for AFLA-maize and FER-maize was clearly demonstrated. This coupled to the large data set used, comprising a 13-year time series, and the good results for the statistical scores applied, together confirmed the robustness of the models developed here.


INTRODUCTION
Mycotoxin contamination of maize is a major concern worldwide (Eskola et al., 2020). The colonization of maize ears by Aspergillus section Flavi and Fusarium spp. can lead to ear rots whose impact on the amount of grain yield is minor or negligible yet their mycotoxin contamination levels are high; therefore, the mains impact of mycotoxin producing fungi in maize regards grain safety and its compliance with the legal limits. Concerning those mycotoxins produced by Aspergillus section Flavi, among the aflatoxins (AFs), aflatoxin B 1 (AFB 1 ) is classified by IARC (International Agency for Research on Cancer) as a class-1A, human carcinogen. Such AFs were first detected in Italy in the early 2000s (Piva et al., 2006;Battilani et al., 2008a), but since 2012 they have spread all over southeastern Europe, presumably aided by warmer and drier conditions during summer attributed to ongoing climate change (Dobolyi et al., 2013;Levic et al., 2013;Battilani et al., 2016), and their incidence and severity can vary markedly among years. Fusarium spp. can produce a wide range of mycotoxins, of which the fumonisins B 1 , B 2 , and B 3 (FBs)-predominantly produced by F. verticillioides-are the key ones reported in maize grain worldwide, thus posing a serious risk of possible human carcinogenicity (IARC, 1993). Other mycotoxins produced by Fusarium genus are the trichothecenes (TCTs) and zearalenone (ZEN), these being prevalent in temperate and wet areas especially in rainy years, optimal conditions for their main producer, F. graminearum (Pietri et al., 2004).
These main mycotoxins threaten the maize supply chain worldwide, in all producing areas; nevertheless, the prevailing mycotoxin and level of contamination depends both on the growing area and year, intended as the meteorlogical conditions occurring during the crop growing season (Logrieco et al., 2021). Support for farmers coming from predictive modeling, using meteorological data as input variables (Battilani, 2016), has been pursued in Europe in the form of two mechanistic models for AFB 1 and FBs predictions: respectively, AFLA-maize (Battilani et al., 2013) and FER-maize (Battilani et al., 2003). They both aim to predict the risk of contamination above current legal limits in force in Europe (European Commission, 2006b, 2007, and they strongly support stakeholders in the maize chain management Battilani, 2016;Palumbo et al., 2020). However, mounting uncertainty of climate conditions and extreme events, often emphasized as issues in climate change, has recently increased the importance of deriving reliable predictions at the farm level (Camardo Leggieri et al., 2020a). Addressing the variability in mycotoxin occurrence among years and geographic areas, even those quite close to each other, in addition to the emerging issue of co-occurring mycotoxins (Camardo Leggieri et al., 2019;Giorni et al., 2019), will require making reliable predictions to support the maize chain management.
Weather variables are the leading factors contributing to mycotoxin occurrence, but the cropping system used is a powerful tool of farmers to mitigate grain contamination. Accordingly, several authors have studied the role of the cropping system and the rationale behind its impact on mycotoxin contamination (Munkvold, 2003;Battilani et al., 2008a;Blandino et al., 2009;Kos et al., 2013;Palumbo et al., 2020). A rationale crop rotation, leaving the land fallow (unsown with maize), is recommended to reduce mycotoxin contamination in maize fields, even if the impact on it cannot be readily demonstrated, especially in intensive maize-growing areas (Guo et al., 2005;Marocco et al., 2008;Munkvold, 2014). A significant impact of the season length of maize hybrids, frequently reported as FAO class (Food and Agriculture Organization classification), upon FBs and AFB 1 contamination was reported . Scientist do not all agree on this statement (Battilani et al., 2008a;Mazzoni et al., 2011), but the number of days elapsed from sowing to harvest was positively related to mycotoxin contamination (Battilani et al., 2008a;Torelli et al., 2012). The sowing date has been confirmed to influence the likelihood and extent of mycotoxin contamination (Jones, 1981;Alma et al., 2005), with late sowing generally associated with a higher content of mycotoxins at harvest (Alma et al., 2005;Blandino et al., 2009;Mazzoni et al., 2011). Nonethless, irrigation has a strong impact as well, particularly upon AFs occurrence (Palumbo et al., 2020). Both the severity of European corn borer (ECB) (Jones, 1981) and the use of insecticide treatments may also significantly affect contamination, especially from FBs (Alma et al., 2005;Saladini et al., 2008;Mazzoni et al., 2011). Harvest time, or rather the kernel moisture at harvest, can also be crucial, notably for AFs contamination of maize (Munkvold, 2003;Battilani et al., 2008a); in fact, AFs production increases significantly from maize physiological maturity, when kernel moisture is lower than 28-30% (Payne et al., 1988;Giorni et al., 2016). Therefore, delays in maize harvest after that stage means giving the fungus time to efficiently increase the contamination. Then, keep kernel moisture below 14% is mandatory in the postharvest stages, from drying to the whole storage period (Danso et al., 2018).
The above research findings have contributed to developing guidelines for mitigating mycotoxin contamination, but a quantitative evaluation of cropping system's impact on mycotoxins remains an unresolved issue. The little work done so far to predict the effect of cropping system on mycotoxin contamination (Battilani et al., 2008b;Camardo Leggieri et al., 2015) is neither complete nor satisfactory; however, it does clarify that cropping factors cannot be considered in isolation and that applying conventional statistical methods is not suitable for the task (Battilani, 2016;Palumbo et al., 2020). Hence, alternative approaches should be explored and possibly used.
Machine learning's emergence, alongside big data technologies and high-performance computing, introduces new opportunities for data-intensive science in precion farming and sustainable agriculture (Liakos et al., 2018). ML is the scientific field in which machines are trained to learn without being strictly programmed (Samuel, 2000), which has three main categories: (1) supervised learning (SL), (2) unsupervised learning (UL), and (3) reinforcement learning (RL). The SL algorithms use a training data set of labeled data to infer a function that is used to map new data. The UL algorithms directly look at the data and learn patterns from them, without human supervision. Last is RL, the ML branch in which entities called "software agents" take action, in a specific context, to optimize a given function. These ML approaches are increasingly applied in different subject areas to solve complex problems, often those with many factors involved, to which agriculture is no exception. In fact, ML is used in a variety of contexts and all the three main categories are now applicable (Liakos et al., 2018;Elavarasan and Vincent, 2020). Recently, Liakos et al. (2018) reviewed the ML approach in agriculture, highlighting that ML models had been applied in the multi-disciplinary agri-technologies domain for crop management (61%), yield prediction (20%), and disease detection (22%), but never accounting specifically for mycotoxins' cooccurrence. In crop yield prediction, which depends on many different factors operating simultaneously, deep neural networks (DNN), a type of artificial neural network (ANN) for SL models, are the most used Nevavuori et al., 2019;Niedbała, 2019). DNNs are also very useful for plant disease identification, which is done via convolutional neural networks (CNN), which is a specific DNN-architecture used for image recognition (Boulent et al., 2019). Mycotoxins are mainly detected via high-performance liquid chromatography (HPLC) and mass spectrometry; however, DNNs coupled with rapid analytical tools, such as the electronic nose or infrared attenuated total reflection spectroscopy, have been recently applied and found to improve the assessment reliability (Evans et al., 2000;Jia et al., 2019;Öner et al., 2019;Camardo Leggieri et al., 2020b). Torelli et al. (2012), in the first example of ML applied to mycotoxins, performed a 2-year study (2007)(2008) that included seven cropping system variables-FAO class, sowing and harvest dates, crop duration, kernel moisture, ECB treatment, and irrigation-as input for an ANN to classify maize samples based on their contamination with FBs. A fair correlation between the predicted and observed contaminated samples was reported (R 2 = 0.67 and R 2 = 0.57, respectively, for the training and validation data sets), so the approach seems promising. Therefore, this study aimed to develop ML models by combining the AFLA-maize and FER-maize predictions, together with cropping system information, as input to improve the mycotoxin risk predictions for AFB 1 and FBs. For this, A 13-year data set was considered and two ML models, one for each type of mycotoxin, were trained and validated (5-fold cross-validation) and their respective performance discussed.

Data Collection
The data used in this study came from several surveys managed in the Emilia Romagna region (northern Italy) during 2004-2018, partially published (Battilani et al., 2013;Camardo Leggieri et al., 2015, 2020b. The protocol for data collection was the same both for the published and unpublished data. Briefly, meteorological data were downloaded from the Emilia Romagna meteorological service, available on request for research applications. They were based on a grid of squares, each 5 km wide, that encompassed the Emilia Romagna territory; all sources (both meteorological stations and radar) are interpolated for each square, to deliver reliable data (Bottarelli and Zinoni, 2002). Hourly data on air temperature (T, • C), relative humidity (RH, %), and rain (R, mm), during the period of January through September, were downloaded. These squares were filtered, to locare those corresponding to the maize field site sampled.
Maize fields sampling was performed at harvest, managed between mid August and September, during the combine machine discharge, according to European Commission Regulation (EU) 401/2006 (European Commission, 2006a). Relevant cropping system data were collected in each georeferenced field, based on a questionnaire filled by farmers, supported by extension services. Empirical information of different site variables were collected: the type of soil (percentages of sand, clay, and silt), maize hybrid FAO class, preceding crop, type of tillage, sowing date, plants per m 2 , silk emission date, harvest date, damage caused by hail or wind and ECB, fertilization type and dose, the number of irrigation intervention with relative volumes of water used, pest and disease control practices, and kernel moisture at harvest. Mycotoxin analysis was performed for all the sampled fields according to Bertuzzi et al. (2012) for the AFs [limit od detection (LOD): 0.05 µg/kg and limit of quantification (LOQ): 0.15 µg/kg], and by following Pietri and Bertuzzi (2012) for the FBs (LOD: 10 µg/kg and LOQ: 30 µg/kg).

Data Processing
AFB 1 and FBs content allowed in maize grain, according to legal limits, were used as a threshold to separate the field samples in two classes: (1) contaminated, consisting of those samples equal or exceeding the respective legal limit; (0) noncontaminated comprising all samples below the legal limit. Thresholds were therefore set to 5 µg/Kg for AFB 1 and 4,000 µg/Kg for FBs (FB 1 + FB 2 ), the legal limits currently set by the European Commission for unprocessed maize destined for human consumption (European Commission, 2006b, 2007. Meteorological data were used as input for the two predictive models, AFLA-maize and FER-maize, and cumulative risk indexes were obtained as the output, AFI for AFB 1 and FK for FBs, for each station and year. DNN models were implemented within the frame of Scikit-Learn (v0.21.3) in the Python module library (Pedregosa et al., 2011).

Input Features
After the exclusion of those variables with many missing data points, eight different variables were considered as input for the ML approach; sowing date and harvest date were grouped on a per week basis. Of these eight, five variables were categorical (maize hybrid FAO class, preceding crop, sowing week, harvest week, ECB damage; Table 1) and three were continuous variables (growing days, days from crop sowing to harvest, calculated variable, kernel moisture at harvest, and mycotoxin cumulative indices AFI and FK, the output of predictive models).
For each continuous variable, its average (µ) and standard deviation (σ) were computed for the data standardization, using Eq. (1). Even if this procedure is entirely facultative in neural network training, it is nonetheless useful for reducing the variance, speeding up the computational process, and improving the model's accuracy (Jin et al., 2015).
To deal with categorical data, the integer encoding procedure was applied. This assigns to a specific category an integer value that ranges from 1 to N, where N is the last category.

Deep Neural Network (DNN) Development
A typical ANN consists of a network of connected computational units called neurons; these units are organized in layers, with input data passing through the network, and an activation function used to produce an output. DNNs are a particular class of ANNs in which, between the input and the output layer, there is an arbitrary number of hidden layers. The fully connected architecture has been adopted, meaning that each neuron in each layer is connected to each neuron in the next adjacent layer.
The development of a DNN model is a two-step process: (i) training and (ii) validation. The final aim of the training is to minimize a given error function by using an optimization algorithm. The training phase ends when the error converges to a pre-determined value, or when it does not decrease for a specific number of cycles, both decided a priori by the user (Camardo Leggieri et al., 2020b). A "batch training" mode was applied in this work. Briefly, during the training phase, training data (i.e., the mycotoxin contamination data in this study), were split into subgroups called batches. Neural network weights are updated when every sample inside a batch passes through the network. The iteration ends when all batches have passed through the network. After each iteration, a penalization term, called the weight decay (L2 regularization term), is introduced into the model to avoid overfitting the model to the data.
The final output is the result of a linear or non-linear activation function. In our study, a non-linear activation function between the input layer and the hidden layers, called Rectified Linear Unit (ReLU, Eq. 2), was applied: where, x represents the weighted sum in a given input to a neuron (Dahl et al., 2013;Zeiler et al., 2013;LeCun et al., 2015). The activation function used between the last hidden layer and the output layer (classification function) took the logistic form (Eq. 3): where, j is relative to the j th output neuron, and i represents the i th input neuron. The numerical result is between 0 and 1, for which 0.5 served as a threshold to discriminate between the positive and negative classes. Different DL models were tested following a grid search procedure, done as described in Camardo Leggieri et al. (2020b).
Briefly, each hyperparameter was tested with a combination of every other hyperparameter. The hyperparameters tested were weight decay, number of hidden layers, number of neurons per hidden layer, and the optimization algorithm. In this work, two algorithms where tested; the first approximates the Broyden-Fletcher-Goldfarb-Shanno algorithm and is called LBFGS (Byrd et al., 1995;Kingma and Ba, 2014), while the second is an optimization of the classical stochastic gradient descent, called Adam (Kingma and Ba, 2014). Matthew's correlation coefficient (MCC) and accuracy were used as metrics to select the best combination of hyperparameters, for both NN models relevant to the pathosystem A. flavus-maize (DNN-A. flavus-maize) and F. verticillioides-maize (DNN-F. verticilloides-maize).

DNN Validation
As in Camardo Leggieri et al. (2020b), both AFB 1 and FBs original datasets were splitted into two "sub-dataset." These four "sub-datasets" (two for AFB 1 and two for FBs) were generated by random sampling, but keeping the proportion of contaminated vs. non-contaminated samples constant. The first "sub-data set" accounted for 75% of the original data set and was used to perform a 5-fold cross-validation (CV). The other, called the "blind set" in this work, accounted for the rest of the original data set (25% of the original data set), and was used for a further validation of the models. The goodness-of-fit of each DNN-A. flavus-maize and DNN-F. verticilloides-maize model was assessed by computing several statistical scores, which were also applied to the blind set: • True positive and true negative rates (TPR, TNR; Kohavi and Provost, 1998); where, TP and TN denote the number of true positives and true negatives, respectively, and likewise FP and FN denote the number of false positives and false negatives.
• Positive predictive value (PPV): the PPV (Eq. 5) index represents the proportion of positives samples identified as true positives. It ranges from 0 to 1 (Kohavi and Provost, 1998) • Receiver operator characteristic (ROC) curve and area under the curve (AUC): the ROC curve and its AUC measure the quality of a binary classifier: the higher the area under the curve, the better the model performs. The AUC value ranges from 0 to 1 (Bradley, 1997). • The MCC (Eq. 6) is used to assess the quality of binary classification (Matthews, 1975). This index takes into consideration TPR, TPN, and both false discoveries (false positives and negatives). The MCC ranges between -1 (complete disagreement between predicted and observed values) and +1 (perfect agreement). The MCC is considered a balanced measure, and it can be used even if the two classes differ in size (Boughorbel et al., 2017).
All the scores were computed using a home-built Python (v3.6.9) script that implemented the equations reported above. The ROC curves and AUCs, the grid search procedure, the MCCs, and the DNN architecture were implemented in the framework of scikitlearn (v0.23.2; Pedregosa et al., 2011). The output indexes of the two mechanistic models (AFI and FK) were used to classify the blind data set, as described in Battilani et al. (2003) for the FBs and in Battilani et al. (2013) for AFB 1 . Finally, the results were compared using the classification obtained by the two DNNs.

RESULTS
A total of 378 and 225 samples were included in the A. flavusmaize and F. verticilloides-maize data sets, respectively ( Table 2). No data were retrieved for FBs before 2009 or during 2012 and 2013. Concerning the AFB 1 , data for it were not retrieved for the years 2012 and 2013. The sample sizes per year were slightly different for the two data sets because of this missing data.
Fields with AFB 1 contamination above 5.0 µg/kg were found, whose incidence and mean values differed across the considered years. The highest amount of AFB 1 was found in 2009, at 494.3 µg/kg. In all years AFB 1 was greater than LOQ, but lower than LOD, except in 2009, 2014, and 2016. Regarding the incidence of fields found positive for AFB 1 , this was highest at 50% in 2017 and the lowest (12.9%) in 2011 ( Table 2). In the whole maize data set, for AFB 1 , the mean incidence of positive samples was 32.9%.
Fields with FBs' contamination above 4,000 µg/kg were found in all the years considered, but this incidence differed across all years. The only year when the FBs was below the LOD was in 2017. The highest amount of FBs (106,053.5 µg/Kg) and the highest incidence of positive samples (89.8%) were both detected in 2014. The lowest incidence of FBs was 6.7%, scored in 2011 ( Table 2). Considering the FER-maize data set as a whole, the overall mean incidence of positive samples, above the legal limit, was 32.2%.
Means and standard deviations computed for AFI and FK, the kernel moisture, and the growing days, are reported in Table 3.
Categorical data were grouped into three or four categories ( Table 1). The FAO class was represented by four categories: 200-300, 400, 500, and 600-700. The preceding crop was represented by three categories: arable crops, small grain, and maize. The sowing and harvest weeks both accounted for four categories. Considering the sowing week, category #1 was assigned to weeks 10-12 of the year, #2 to weeks 13-14, #3 to weeks 15-16 and #4 to week ≥ 17. For the harvest week, those weeks of the year from 32 to 35 were designated category #1, and likewise 36-37 to # 2, 38-39 to #3, with all weeks > 40 assigned to #4. The damage caused by ECB was divided into three categories: no damage and small damage were grouped into category #1, medium damage was assigned to category #2, and severe damage was assigned to category #3 (Camardo . Standardized continuous data were joined to categorical data, to form the neural network's input vector; thus, the final input array was formed by five encoded and three continuous variables (Tables 1 and 3), respectively.
The generation of the "sub-datasets" included a total 283 (75%) and 95 (25%) samples in the CV and blind set respectively

Deep Neural Network (DNN) Training
The two DNNs were trained to be able to predict the content of AFB 1 and FBs, respectively. The NN-A. flavus-maize model was implemented with one hidden layer consisting of 80 neurons, a ReLU activation function, and an L2 regularization term of 0.0001. The parameters of that model were updated using the Adam algorithm (Kingma and Ba, 2014). By contrast, the NN-F. verticilloides-maize model's sole hidden layer had 50 neurons, a ReLU activation function, and an L2 regularization term of 0.1. Its parameters were updated using the LBFGS algorithm (Byrd et al., 1995; Table 4). The two developed NN models were validated using both the 5-fold cross-validation and a blind data set.

DNN Validation
Cross-validation ROC curves and their relative AUCs were computed for the two NN models, to assess the quality of the two classifiers (Figures 1A,B). The 5-fold-cross-validation for NN-A. flavus-maize achieved an accuracy of 66.56 ± 3.381 (mean ± SE), and even higher at 78.94 for the blind data set. An MCC of 0.10 ± 0.157 and 0.49 were achieved by the crossvalidation and blind data set, respectively. The AUC and the TPR averaged 0.58 ± 0.063 and 0.08 ± 0.073 during the model's crossvalidation. The model scored an AUC of 0.64 and a TPR of 0.42 when tested against the blind data set ( Table 5).
The NN-F. verticilloides-maize model's 5-fold cross-validation attained an accuracy of 69.63 ± 10.892; the accuracy was 79.31% using the independent data set. A TPR of 0.53 ± 0.118 and 0.65 were, respectively, achieved by cross-validation and blind data set, respectively. Moreover, the model had an AUC of Parameters update algorithm

Adam LBFGS
The number of input neurons is defined as the sum of the categorical and continuous data. The hyperparameter values were selected by following a grid searching procedure.
0.72 ± 0.103 and an MCC of 0.35 ± 0.229 during the crossvalidation phase, with corresponding values of 0.75 and 0.56 when tested against the unseen data set (Table 4). Finally, to check whether the new approach represented a major step forward in the prediction of mycotoxin contamination in maize, our two DNN-models were compared with two counterpart mechanistic models, both run with the whole available data set. The resulting confusion matrix (Table 6) shows the performances of the NN-A. flavus-maize and NN-F. verticilloides-maize models and those of the two mechanistic models vis-à-vis the blind data set. The NN-A. flavusmaize model correctly estimated about 78% of samples (14% true positives, 64% true negatives). The wrong classification accounted for 19% of them being underestimated and 3% overestimated. The NN-FER-maize model correctly classified approximately 75% of the data set (25% true positives, 50% true negatives), with underestimations and overestimations amounting to 15 and 11%, respectively. In stark contrast, the AFLA-maize model correctly predicted just ∼53% of samples (11% true positives, 42% true negatives). Further, a wrong classification accounted for more underestimations (22%) and overestimations (25%). Similarly, the FER-maize model correctly classified only ∼52% of samples (31% true positives, 20% true negatives), but its wrong classifications included fewer (7%) underestimated cases being more prone to overestimations (41%).

DISCUSSION
Maize is exposed to mycotoxins, which threaten human and animal health, and represent the major non-tariff trade barrier for agricultural products, negatively affecting the income of smallholder farmers and disrupting regional and international trade (Palumbo et al., 2020;Logrieco et al., 2021). Timely identification of contaminated lots is not a trivial challenge since mycotoxin contamination relies on several factors, including meteorology and how farmers manage the crop during the season and in the postharvest stages of storage and distribution (Munkvold, 2014;Logrieco et al., 2021). Different methodologies for the rapid detection of mycotoxin contamination are currently available (Öner et al., 2019;Camardo Leggieri et al., 2020b), but since they are applied at harvest or postharvest stages, they offer no support for taking preventive action and for optimizing lot use and management. On the contrary, farmers can benefit from model predictions of the risk of mycotoxin occurrence above the legal limit, when delivered before or during the cropping season, in the form of risk maps or risk indexes. Therefore, predictive modeling has garnered mounting interest over the last two decades (Battilani, 2016). Predictions refer to maize at harvest and it is assumed that the postharvest management guarantee a rapid grain drying to humidity ≤14%, kept stable during storage, to avoid fungal activity and further mycotoxin production.
Meteorological factors jointly determine whether fungi can grow and produce toxins, while the site's cropping system modulates the amount of contamination that ensues . The former are the driving  The deep neural networks (DNNs) considered were evaluated using a 5-fold crossvalidation and a blind data set. Cross-validation (CV) results were reported as an average value and related standard deviation (Avg ± Stdev). ACC, accuracy; MCC, Matthew's correlation coefficient; AUC, area under the curve; TPR, true positive rate; TNR, true negative rate; PPV, positive predictive value; n.c., not computed.
variables for predictive modeling, whereas the latter are rarely included, especially in mechanistic models. This omission is starting to gravely limit the reliability of predictions; in fact, during the last two decades, the typical cropping system has changed significantly due to the knowledge transfer from scientists to farmers; farmers are now following the guidelines to optimize crop management and mitigate mycotoxin contamination, with good results so far in term of a reduced mycotoxin occurrence. Both the meteorological data and the cropping system data have been used before as model inputs, to predict the content of mycotoxin in maize at the time of its harvest (Battilani et al., 2003(Battilani et al., , 2008aBertuzzi et al., 2014;Camardo Leggieri et al., 2015), yet they were used independently and only supported by basic statistical approaches. The aim of this study was to evaluate how combining cropping system information with mechanistic predictive models could support the sought-after improvement in prediction performance.
Here an ML approach was developed using the AFLA-maize and FER-maize outputs (mycotoxin risk indexes) combined with cropping system information-this being known to significantly influence mycotoxin contamination in maize according to other studies (Palumbo et al., 2020;Logrieco et al., 2021)-as input variables. Other crop-related variables should have been included, like fertilization, irrigation, and pest control (Mazzoni et al., 2011;Munkvold, 2014); however, we excluded them because this data was largely unavailable to us. Moreover, the geolocation of maize fields was excluded as an input variable in our modeling; actually, even when the field location is known to be relevant (Torelli et al., 2012;Camardo Leggieri et al., 2015), the idea of this work was to obtain models applicable at a global level, without geographical constraints. When combining weather data and cropping system no information is lost, even when the maize fields' geolocation is excluded.
The predictions of NN-A. flavus-maize and NN-F. verticilloides-maize, the two neural network models developed in this study, were capable of an accuracy approaching ∼78% for AFB 1 and ∼79% for FBs, with a good correlation between Both AFLA-maize and FER-maize, the mechanistic models which served as the starting point of our investigation, achieved accuracies one-third lower, of about 50%, and their respective MCC was very close to 0; this indicates that their predictions were comparable to random one, when based on the same data set (i.e., the blind data set) used for NN model evaluation. It is, therefore, evident the proposed ML approach significantly improved the prediction of mycotoxins' content across the studied maize fields, making the successful use of this tool to detect maize grain not compliant with the current legal limit in Europe now more realistic and feasible to implement. Nevertheless, the NN-F. verticillodes-maize model performed better than NN-A. flavus-maize; apparently, it is easier to predict levels of FBs than AFB 1 as contaminants. The reason for this is not entirely clear, but the very low limit fixed by the legislation for AFB 1 could surely play a role. Further, irrigation is known to be very relevant for AFB 1 contamination, and grain humidity at harvest too, but they were excluded as input variables because of many missing data and this has surely a considerable impact on the predictive capacity of the model. The NN-models were developed based on a large data set, one that included data collected from over 12 different years, with 378 samples in the AFB 1 data set and 225 samples in the FBs data set. The large data sets used, and the range of years considered, including the period of significant changes in the cropping system, strongly support the robustness of NN-A. flavus-maize and NN-F. verticilloidesmaize models and their promising utility as a tool to support farmers in their decision-making. Future applications in other pathosystems is also foreseeable, as previously done for AFLA-maize (Kaminiaris et al., 2020).

CONCLUSION
To conclude, despite omitting some relevant cropping system variables, a substantial improvement at correctly predicting maize fields contaminated with mycotoxins above their legal limits was gained. Further improvements should be obtained by optimizing the data collection. Solving the missing data problem might be an easy task in the future, once scientists succeed in convincing farmers of the crucial role they can play in such data collection and also of the added value of predictive models in mycotoxin management. This will be a matter of building with the maize chain stakeholders a knowledge exchange approach and make them more involved compared to what actually happened in the past. Another aspect of improvement could be gleaned through our modeling approach, that of emerging issues related to mycotoxins. The main example of this concerns the recent report of co-occurring A. flavus and F. verticilloides in maize ears due to climate change effects, resulting in their complex interaction, with their dominance alternating during the growing season (Camardo Leggieri et al., 2019Leggieri et al., , 2020aGiorni et al., 2019). Looking ahead, we anticipate that elucidating the impact of these interactions between coexisting fungi upon mycotoxin production in maize will become crucial. Data collection to develop a joint model for the prediction of AFB 1 and FBs, including the impact of fungi co-occurrence, is ongoing and in the next future it is expected to contribute to a step up in mycotoxin prediction, possibly joining also NN-A. flavus-maize and NN-F. verticilloides-maize in a NN-mycotoxmaize predictive model.
The current work represents a notable step forward in modeling and predicting mycotoxins in crops. We retrieved evidence that ML can effectively combine cropping system data and meteorological data, thereby improving the accuracy and robustness of predictions. Big data is a relatively new concept for agriculture and plant disease research and management, but massive volumes of data with several components that interact within the pathosystem can also be captured in this context, and their elaboration can enhance the decision-making process (Wolfert et al., 2017). Applying machine learning to farm management systems is quickly evolving into a real artificial intelligence (AI) system, providing richer recommendations and insight for subsequent decisions and timely actions (Liakos et al., 2018). Further research that aims to integrate automated data recording, mycotoxin analysis, ML implementation, and decision support systems will provide practical tools in line with so-called "knowledge-based agriculture." This should move us closer toward sustainable agriculture and smart farming that also improves food safety and quality.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
PB, MC, and MM contributed to the conception and design of the study and wrote sections of the manuscript. MC organized the database. MM performed the statistical analysis. MC and MM wrote the first draft of the manuscript. All authors contributed to manuscript revision, read it, and approved the final submitted version.

FUNDING
This study was partially supported by the regional government, Emilia Romagna region, "SERVICE -SistEmi infoRmatiVi rIschio miCotossinE" (IT systems for mycotoxin risk) No. 5149128.

ACKNOWLEDGMENTS
We are grateful to the Emilia Romagna meteorological web service for providing meteorological data, to CRPV for the project coordination and to all the farmers and technicians who contributed to maize samples collection.