We conducted spatial–temporal monitoring to investigate the occurrence of A. catenella. The monitoring site is located on the southeastern coast near Namhae, Geoje Island, and Busan, South Korea (Figure 1), opened toward off sea, causing fresh and oceanic water intrusion (Kang et al., 2012). In particular, the eastern coast of Geoje Island has frequently occurring PSP initiation (National Fisheries Research and Development Institute, 2020). PSP outbreaks provoked by A. catenella bloom in our study area have been reported (Chang et al., 1987; Kim et al., 2015). Water surface elevation, wind, and wind velocity data were measured by the Korea Ocean Observing and Forecasting System and used to set up the model. Two observations at stations 1 and 2 were used to calibrate the water level, temperature, and salinity. Another three observations at stations 3, 4, and 5 were used to calibrate NH4-N and PO4-P (Figure 1).

One-thousand one-hundred and seventy-five samples were collected from the water surface at the sampling sites, which had water depths of 12 to 59 m. The monitoring period was from January 2017 to December 2019. Water samples were acquired with a Van Dorn bottle and fixed with Lugol solution (final conc., 2%) from 9:00 AM to 4:00 PM. The fixed seawater was concentrated to 5–50 mL aliquots by overnight sedimentation. Alexandrium catenella cells were enumerated using a Sedgwick–Rafter counting slide at a 200× magnification with a light microscope (Zeiss Axioscope 2). To identify the species based on morphological and molecular analyses, the cells were processed as described by Kim et al. (2020). Cysts of A. catenella were isolated from sediment samples collected using a core sampler and incubated at bottom temperatures on the sampling dates. The germination ratio of the cysts was estimated by the percentage of cysts germinating within one month. Growth experiments of A. catenella were conducted under several temperature and salinity conditions, reflecting the sampling site environments. The measurements of germination ratio and growth rate were described in detail by Kim et al. (2020).

The input data in this study consisted of two parts: (1) ocean physical data (e.g., water velocity, water temperature, water elevation, and retention time) and chemical data (e.g., salinity, PO₄-P, and NH₄-N) and (2) ocean biological data [e.g., germination ratio, growth rate, and operational taxonomic unit (OTU)]. This study applied the environmental fluid dynamics code (EFDC) model, adopted for simulating ocean and coastal waters (Dai et al., 2011; Du and Shen, 2016) to generate the ocean physical and chemical data. The EFDC model is a fluid simulation model that includes three-dimensional flow and biochemical transport in the ocean, estuaries, and lakes. EFDC can solve free surface, vertical hydrostatic, and turbulent equations for fluids with different densities. The governing equation was derived using the vertical hydrostatic boundary with turbulent equations and consists of the momentum [Eq. (S1, S2)], vertical hydrostatic pressure (Eq. S3), and continuity equations [Eq. (S4, S5)] (Jeong et al., 2010). The ocean physical and chemical data included the temporal and spatial distributions of water velocity, water temperature, water elevation, retention time, salinity, PO₄-P, and NH₄-N. These data have been verified as variables that influence the life cycle of A. catenella (Itakura and Yamaguchi, 2001; Kim and Yoo, 2007; Armi et al., 2011; Kim et al., 2020). The ocean physical and chemical distributions were calculated by the EFDC model. Ocean biological data included the temporal and spatial distributions of the germination ratio of A. catenella cysts, growth rate of vegetative cells of the species, and OTU of bacteria. The cyst germination and growth rates are critical elements in recurrent outbreaks of A. catenella blooms in situ, although dormant populations are significantly affected by environmental variables. Hence, we adopted these rates as the input data to simulate A. catenella. The data pertaining to cyst germination and growth rate were obtained by Kim et al. (2020). The microbial community data in this area, including the OTU data, were extracted from the previous study (Cui et al., 2020). operational taxonomic units were analyzed at a distance of 0.01 using the mothur v1.39.3 pipeline with SILVA database, release 132 (Schloss et al., 2009; Kozich et al., 2013). Three representative OTUs (OTU1, OTU2, and OTU3), which were identified as A. catenella-related OTUs based on ecological network analysis in the previous study (Cui et al., 2020), were selected for further analysis in this study. Taxonomically, OTU1, OTU2, and OTU3 were assigned to genera Fluviicola (family Crocinitomicaceae), Ascidiaceihabitans (family Rhodobacteraceae), Candidatus Actinomarina (family Candidatus Actinomarinaceae), respectively. Supplementary Figure 1 presents the process used to generate the temporal and spatial distribution of biological data. A linear model was adopted to generate the variation of biological data by changing the environmental variables. The linear regressions of the germination ratio and growth rate were calculated using temperature and salinity from in vivo experiments (Supplementary Figure 1A), whereas that of the OTU distribution was generated using the temperature, salinity, PO₄-P, and NH₄-N values from monitoring (Supplementary Figure 1B). A linear regression model explains the relationship between one response variable and multiple explanatory variables (Montgomery et al., 2021); therefore, these linear models used the simulated water temperature, salinity, PO₄-P, and NH₄-N from EFDC to determine the temporal and spatial distribution of biological data (Supplementary Figure 1C). The simulation period was from January 2017 to December 2019. These ocean physical, chemical, and biological data were validated using observational data. The EFDC grid was a Cartesian grid with a cell size of 200 m × 400 m. Observations of wind direction, velocity and water surface elevation were used to set up the EFDC model. More details of the experiments on germination ratio, growth rate, and OTU are presented in the Supplementary Information.

Figure 2 shows the two-step process of generating input data to apply CNN for simulating the bloom of A. catenella: (1) extracting the spatial–temporal distribution of the physical, chemical, and biological information at the study site (Figures 2A,B) and (2) converting the input data to two-dimensional data (Figure 2C). These physical, chemical, and biological data have three dimensions, m × m × (n + 1), and the simulation point of A. catenella is located in the center. Here, m is the size (i.e., height × width) of the input window that includes the simulation point, and n is the lookback, which indicates the number of time steps included based on the current simulation time. For example, when m = 3 and n = 5, the data size is 3 × 3 × (5 + 1). Because the CNN approach was developed focusing on two-dimensional data (length × height) or RGB data (length × height × 3), we converted these inputs into two-dimensional data before applying the CNN (Shin et al., 2016; Koundinya et al., 2018).

The simulation layout of A. catenella is presented in Figure 3. The simulation comprises of the following three steps: (1) optimizing the input and CNN structures, (2) simulating the bloom of A. catenella based on the optimal input and CNN structures, and (3) identifying A. catenella bloom factors based on Grad-CAM. Because these CNN structures (e.g., Resnet, GoogLeNet, and Inception) were developed using data of size 299 × 299, the input data had to be converted (Szegedy et al., 2016; Akiba et al., 2017). Subsequently, the converted data are fed to the classification and regression CNN models; the classification CNN model decided the initiation of A. catenella, while the regression CNN model generated the density after A. catenella was initiated. The initiation and the density indicated the occurrence and the number of A. catenella cells, respectively. The use of two CNNs with different roles can prevent bias in model training, as most of our monitoring data were zero, indicating that A. catenella did not occur. In addition, we optimized the input window (m), lookback (n), and CNN structures (Figure 3A) to generate the spatial–temporal distribution of A. catenella (Figure 3B). Using the optimal parameter values and CNN structures, we analyzed the performance of A. catenella forecasting with increasing forecast lead times (days) of up to seven days. Model training was performed using Intel ^® Xeon CPU E-52687W v4 @ 3.00 GHz, 128 GB RAM, and NVIDIA GTX 1080 Ti. CNN was implemented with the machine and deep learning toolboxes in MATLAB. Accuracy was used for evaluating the classification CNN model, while RMSE and R² were used for the regression CNN model. Relevant explanations can be found in the Supplementary Information.

Convolutional Neural Network (CNN) is a popular deep learning model that extracts data features using convolving filters (Deng et al., 2009). A typical CNN architecture consists of convolutional, pooling, ReLU, batch normalization, concatenation, normalized, and fully connected layers (LeCun et al., 2015). Each layer has a specialized role in the architecture: convolutional and pooling layers are used for feature extraction. ReLU and normalization layers are used for a linear and normalization calculation, respectively (LeCun et al., 2015). These layers can be combined to enhance model performance (Szegedy et al., 2016; Khan et al., 2019). Details concerning the CNN layers can be found in the Supplementary Materials document. Supplementary Figure 2 shows the CNN structures employed in this study. GoogleNet and Inception v3 models adopt parallel CNN layers (Szegedy et al., 2016). ResNet 50 and ResNet 101 have skip connections, i.e., the information of the output layer is transferred into the next layer, and into the earlier layer (Szegedy et al., 2016). This structure can increase model performance by reducing overfitting and vanishing gradient problems (He et al., 2016a). In our study, the softmax and mean squared error (MSE) (see Supplementary Information for more details) were used as the loss function for classification and regression CNN, respectively. The loss functions were applied for calculating the error between simulation and observation during model training. For model training, the CNN used hyperparameters, including epoch number, batch size, and learning rate (Loussaief and Abdelkrim, 2018). Epoch is the number of times the learning worked in the entire dataset, whereas batch size is the number of samples used for training (Robert, 2014). The learning rate controls the step size at each iteration to minimize the loss function (Robert, 2014). The CNN comprised 500 epochs with a mini-batch size of 32; the applied learning rates equaled 0.001 and 0.0001 for the first 200 and remaining 300 epochs, respectively. Each epoch generated a corresponding model, and the final model was selected to produce the lowest validation accuracy and MSE. Our study adopted random sampling to divide A. catenella observations into training and validation sets. A uniform distribution was used for the random sampling. Previous studies have also used random sampling with a uniform distribution to divide the data into training and validation sets (Brion et al., 2002; Caruana and Niculescu-Mizil, 2006).

The deep-learning model used in this study applied Grad-CAM for identifying factors contributing to A. catenella bloom (Figure 3C). The Grad-CAM localization map describes the simulation results by highlighting the important regions (Selvaraju et al., 2017). The model interpretability technique is proposed as a strategy that enables the input-based understanding of the results obtained because neural networks are incapable of explaining model results (Selvaraju et al., 2017). Several prior studies have verified the use of Grad-CAM in the visualization of model features (Selvaraju et al., 2017; Chen et al., 2020). This method is based on class activation mapping (CAM) that can extract the significant features by emphasizing the input data region. However, the use of CAM is restricted to CNN structures that comprise a global average pooling layer (Selvaraju et al., 2016). Grad-CAM overcomes this limitation using gradient information from the final convolutional layer for visualizing the important input data regions (Selvaraju et al., 2017). A detailed description of Grad-CAM is presented by Selvaraju et al. (2017).

We compared the simulated and observed water elevation, salinity, water temperature, NH₄-N, and PO₄-P results. The physical and chemical simulations are illustrated in Supplementary Figures 3, 4. The coefficient of determination (R²) of water temperature and elevation was above 0.90 at stations 1 and 2, while salinity was 0.26 and 0.65 at stations 1 and 2, respectively. The average root-mean-squared errors (RMSEs) of water temperature, elevation, and salinity were 0.82°C, 0.07 m, and 1.12, respectively. Compared to water temperature and elevation simulations, the salinity simulation had a larger error than observation. The average coefficient of determination (R²) of NH₄-N and PO₄-P were 0.67 and 0.30, and the average RMSE values were 0.09 and 0.01 mg L^–1, respectively. Although both NH₄-N and PO₄-P simulations follow the observation trend, their values were underestimated at the peak point. The linear regression equations of growth rate, germination rate, and OTU are presented in Supplementary Figure 5. The growth and germination rates were more strongly influenced by temperature and salinity, respectively, whereas the OTUs related to A. catenella were affected by salinity and PO₄-P.

Figures 4(A.1-4,B.1-4) shows the performance of classification and regression CNN, respectively, with respect to input window size (m) and lookback (n). Here, m denotes the input window size (e.g., height × width), including the simulation point that was located at the center of the input window, and n is the number of time steps. For example, if the input window size is three and lookback is five, the model considers a spatial distribution with 3 × 3 grid cells and temporal information from the previous five days to the current simulation time. In the classification model, except GoogleNet, there existed multiple optimal designs (m, n) in each structure; GoogleNet had an optimal design structure of (3, 30). The model performance of other structures deteriorated when m was above five and n was below fifteen (days). In the regression model, Resnet 50 and Resnet 101 had a similar optimal design of (1, 29), while GoogleNet and Inception v3 required different (m, n) designs; the optimal design of GoogleNet was (1, 2) and that of Inception v3 was (1, 27).

The model performance with optimal input design is summarized in Table 1. Among them, GoogleNet and Resnet 101 offer the best classification performance with an accuracy of 96.83% and RMSE of 1.20 [log(cells L^–1)], respectively. The model performance of Resnet 50 was similar to Resnet 101 with accuracy and RMSE of 96.29% and 1.29 [log(cells L^–1)], respectively. Inception v3 (accuracy of 95.76%) and GoogleNet [RMSE of 1.66 log(cells L^–1)] presented the worst performance in classification (with 95.76% accuracy) and regression models, respectively.

The results of the regression CNN model are plotted against observed A. catenella in Figures 4C.1-4. In Resnet 101, the simulated A. catenella showed good agreement with observations. GoogleNet showed that the simulated A. catenella was overestimated in the low-density cells and underestimated in the high-density cells. The temporal-spatial distribution of A. catenella is presented with observation using GoogleNet and ResNet 101 because these structures showed the best performance in classification and regression CNN (Figure 5A). The simulated distribution substantially followed actual A. catenella blooms. On December 17, 2016, and May 23, 2017, most areas did not provoke A. catenella bloom in both simulated and observed distribution, indicating that our model can simulate this phenomenon. On March 27, 2018, A. catenella blooms were observed in the coastal water. During this period, the simulated distribution of blooms can describe the actual spatial features; the eastern coast presented a relatively higher density than the western one. The data on March 28 and April 25 of 2017 shows increasing A. catenella spring bloom in the study area. The simulated distribution in both periods was in line with the actual distribution; the model generated high-density cells in the eastern coast and non-bloom near Geoje Island of the west coast. On August 19, 2019, the model determined the non-bloom and substantially low density on the western coast. Mismatched results of the spatial distributions are revealed in three spaces without observed data; the channel connected to the northern enclosed bay on May 23, 2017, the western off sea on March 28, 2017, and the enclosed bay on August 19, 2019. Figure 5B shows the performance of A. catenella forecasting with various lead times. All forecast results were found to be worse than those of nowcasting. The average accuracy of the classification model was 95.85% for a lead of up to five days, decreasing sharply thereafter. The average RMSE of the regression model was 1.36 [log(cells L^–1)] until five forecast lead (days) and increased sharply thereafter.

Gradient-weighted class activation mapping (Grad-CAM) shows the feature maps of classification and regression CNN models with GoogleNet and Resnet 101 structures, respectively (Figure 6). These maps were generated depending on the outbreak of A. catenella blooms (e.g., bloom and non-bloom) and density level (e.g., 5–25^th percentile, 25–50^th percentile, 50–75^th percentile, and 75–95^th percentile). The regions with high values are regarded as important features in the map. In the classification model, the important features are affiliated with salinity, temperature, water elevation, latitude-velocity of water, and NH₄-N from 20 to 28 days of lookback when the bloom is not provoked (Figure 6(A.1)). Among these variables, temperature and salinity are the most influenced variables, as evident from the high feature map values. In contrast, the important bloom features highlighted the variables from 3 to 12 days (Figure 6(A.2)). In the regression CNN model, the 5–25^th percentile density show salinity, temperature, water elevation, and NH₄-N from 5 to 30 days of lookback as important variables (Figure 6(B.1)). The important variables in the 25–50^th percentile were similar to the 5–25^th percentile density (Figure 6(B.2)). In the 50–75^th percentile, the lookback from 6 to 27 days is highlighted in the model, while that in the 75–95^th percentiles ranged from 4 to 22 days (Figures 6(B.3), (B.4)) respectively).

The physical and chemical simulation results showed good agreement with the observation results (Supplementary Figures 3, 4, respectively). However, the salinity simulation was worse than that of water temperature and elevation. Previous studies have also shown that salinity cannot improve model accuracy (Hjøllo et al., 2009; Martyr-Koller et al., 2017) because it is vulnerable to external sources (e.g., rainfall and freshwater), increasing simulation uncertainty (Arfib and Charlier, 2016). The simulated water temperature at station 1 was overestimated during winter (December to February), while station 2 provided reliable outcomes in this season. This is because the model was limited to tracking water temperature at two different points if the difference between them was more than 5 °C. Moreover, the NH₄-N simulation followed the observed trends, whereas the PO₄-P simulation was less accurate. Specifically, the simulated NH₄-N was underestimated at the peak point of observation and the simulated PO₄-P was only able to follow the variation in observation. This demonstrated that the simulation of PO₄-P and NH₄-N could not improve model accuracy because these nutrients were observed in low concentrations, making the model sensitive to external sources. Eilola et al. (2009) and Feng et al. (2015) demonstrated that the simulated PO₄-P was underestimated when compared with the observed values, and the simulated nitrogen encountered limitations when following the peak concentration. Additionally, the validation of these simulations is still limited by the lack of observed data. In further research, we will collect more data from additional sites.

In general, the optimal inputs of both the models yielded a small window size below five and long lookback of more than 19 days (Figures 4A,B). Therefore, spatial information below 5 5 grids was sufficient for simulating A. catenella, and temporal information from past 30 days was adequate for improving the simulation performance. The suitable range of grid sizes for the spatial distribution reflected the behavioral characteristics of the species. For the models with window sizes of seven and nine, the simulation performances decreased sharply, indicating that superfluous information might deteriorate the model accuracy. Previous studies have demonstrated that needless data could restrict the effectiveness of model training (Chulkov et al., 2019; Cova and Pais, 2019; Xu et al., 2019). In contrast, in the classification model, the optimal size (m, n) of the input varied depending on the structure based on whether the CNN structures adopted a skip connection or a parallel CNN layer (He et al., 2016b; Szegedy et al., 2016). This demonstrated that the model performance was affected by the structure and the properties of the input data. In addition, a simpler structure, such as GoogLeNet, was appropriate for classifying both bloom and non-bloom (Figure 4(A.3)). In the regression model, Resnet 50 and 101 had similar optimal input properties, whereas GoogleNet and Inception v3 showed lower performances than Resnet. Among the structures, Resnet 101 was the best regression model. This demonstrates that skip connection is appropriate for estimating the density of A. catenella as it can directly connect information from a previous layer to the next layer using skipped layers (Khan et al., 2019). On the contrary, GoogleNet and Inception v3 applied parallel CNN layers having 1 × 1, 3 × 3, and 5 × 5 filter sizes. Hence, the optimization process is important to improve how the input data and model structures affect model training.

ResNet 101 offered the best performance in terms of simulating the cell density of A. catenella that followed the observation distribution, whereas GoogLeNet showed limited performance in simulating A. catenella in low and high densities (Figure 4C). Specifically, the simulated A. catenella from GoogLeNet was overestimated in the case of low density and underestimated in the case of high density, indicating that this model was trained with a narrower range of cell density than that of the other models. The performance discrepancy between ResNet 101 and GoogLeNet may be attributed to the number of layers; ResNet 101 had more layers than GoogLeNet. This demonstrates that the number of layers could affect a model’s performance. According to previous studies concerning CNNs (Khan et al., 2019; Baek et al., 2020), the performance of deep learning models can be influenced by the number of layers. However, a higher number of layers can increase model complexity, leading to the overfitting problem, i.e., the simulated results will correspond too closely to the training data, thereby failing to fit the validation data and predict future simulated results in a reliable manner (Cortes et al., 2017). In contrast, GoogLeNet demonstrated the attainment of a superior validation accuracy (96.83%) as compared to other structures. This confirms that a complex structure is not necessary to identify bloom occurrence and that sufficiently accurate results can be obtained via cell estimation. The simulated distribution substantially followed the feature of actual A. catenella blooms (Figure 5A). Depending on bloom and non-bloom, the model generated similar spatial patterns of density with observation: on December 17, 2016, May 23, 2017, and August 19, 2019, there were no bloom outbreaks in both simulated and observed distributions; on March 27, 2018, the A. catenella blooms were observed in both the distributions; on March 28, 2017, and April 25, 2017, both distributions showed the bloom and non-bloom area in the coast. However, on March 28, 2017, the simulated distribution was determined as blooms in most areas, while there was no bloom of the observed distribution in the western sea. The high uncertainty in population dynamics could be due to in situ environmental and biological variables heterogeneously influencing the initiation, development, and ultimate demise of HABs (McGillicuddy et al., 2015; Brandenburg et al., 2017). The model generated a notable density in the specific spaces without observed data; the simulation of May 23, 2017, and August 19, 2019, presented a higher density in the channel connected to the northern enclosed bay and the enclosed bay, respectively. The mismatch results between simulation and observation also increase the model uncertainty by limiting the quantification of the exact spatial distribution of HABs. Remote sensing can solve this by improving model performance with the spatial distribution of physical, chemical, biological, and atmospheric factors (Shen et al., 2020). Forecast simulation showed lower performance than nowcasting. The performance of the classification and regression model decreased with increasing lead time and deteriorated sharply at a specific date (Figure 5B). These simulation results showed that the model performance would degrade with increasing forecast lead time (days) because the input data might correspond poorly with the future ocean physical and biological processes affecting A. catenella. Prior studies have reported a decline in model performance with increasing lead time. Chattopadhyay et al. (2020) reported a decrease in the model performance from 73 to 47% while predicting a cold-spell class as the lead time changed from 1 to 5 days. Pyo et al. (2020) observed the validation accuracy to decrease with increasing lead time when simulating Microcystis—a causative algal taxon of freshwater HAB. As reported by Miao et al. (2019), an increase in the forecast lead time causes an increase in model uncertainty and imperfect representation of the extracted features. This deteriorated model training and reduced the forecast accuracy. Therefore, this study demonstrated the robust short-term A. catenella forecasting ability of the DL model.

Temperature and salinity are factors with the greatest influence on the classification and regression models, as confirmed by the high values observed in the feature map (Figures 6A,B). Temperature and salinity affect cyst germination and A. catenella growth, thereby causing the proliferation of A. catenella cells (Itakura and Yamaguchi, 2001; Nagai et al., 2004). Ichimi et al. (2001) reported that temperature plays a critical role in cyst germination and algal bloom. Parkhill and Cembella (1999) demonstrated that salinity influences A. catenella growth. NH₄-N and water elevation were also highlighted. A. catenella uses a nitrogen source for growth and prefers ammonium uptake (Siu et al., 1997). Collos et al. (2006) observed that the growth of A. catenella is limited by ammonium uptake and accumulation. Moreover, an increase in water elevation tends to weaken and accelerate the advection and dispersion of algal blooms, respectively (Wu and Kong, 2009). Additionally, Giddings et al. (2014) demonstrated that the physical factors of oceans influence HAB development. The regression CNN model identifies fewer important features that increase the density of A. catenella. In particular, factors, such as the salinity, temperature, and NH₄-N concentration over a 5–30-day lookback period were highlighted in the 5–25^th percentile while higher densities were simulated considering variables over 6–22-day lookback period. The results can be related to the initiation for bloom development of a given inoculum size at a suitable time. The Grad-CAM result reveals that the proposed model generates cell density by changing the highlighted variables. This is because the variables of influence differ depending on the life cycle A. catenella. The water temperature and salinity initiate the development of A. catenella while the nutrient and retention time accelerate its growth (Itakura and Yamaguchi, 2001; Nagai et al., 2004; Armi et al., 2011). Previous studies have reported the transitions in the highlighted region as change input. Panwar et al. (2020) adopted Grad-CAM to identify important areas of COVID-19 detection and demonstrated the change in the highlighted area as an input image. Cheng et al. (2019) extracted hip-fracture features using Grad-CAM and X-ray images. This study estimates the spatiotemporal distribution of A. catenella using CNN and determines the cause of bloom using Grad-CAM. Most previously reported algal simulations have focused on inland water species, such as cyanobacteria, green algae, and diatom. However, only few studies have addressed HAB modeling in seawater and only analyzed limited species (e.g., A. fundyense). HAB modeling in marine ecosystems should consider several environmental factors, such as the temperature, salinity, and velocity of water as well as symbiosis among different species. The paucity of relevant research in this domain could be attributed to the need for the use of complex modeling techniques and availability of limited observations (Ralston and Moore, 2020). This study overcomes this limitation via successful application of deep learning for simulating A. catenella HAB using an ocean model. Further research can be performed considering other species (e.g., Cochlodinium and Karenia) with additional monitoring. The findings of this research are expected to improve the applicability, expendability, and accuracy of HAB modeling using the deep learning models. Therefore, the proposed approach can be considered useful in establishing HAB management in marine environments.